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Abstract  of  Dissertation  Presented  to  the  Graduate  School 
of  the  University  of  Florida  in  Partial  Fulfillment  of  the 
Requirements  for  the  Degree  of  Doctor  of  Philosophy 

FINITE  ELEMENT  ANALYSIS 
OF  COMPOSITE  PLATE  IMPACTED 
BY  A PROJECTILE 

BY 

BRENT  RAY  PETERSEN 
August,  1985 

Chairman:  L.E.  Malvern 

Major  Department:  Engineering  Sciences 

A dynamic  finite  element  computer  program  was  written 
to  determine  the  displacements  and  stresses  in  a layered 
fiber-reinforced  composite  material  plate  impacted  by  a 
projectile.  Yang,  Norris,  and  Stavsky  plate  theory  which 
includes  shear  deformation  and  rotary  inertia  was  used. 

The  variational  form  of  the  plate  equations  of  motion 
was  derived  from  Hamilton's  principle  and  was  used  to 
derive  the  finite  element  plate  equations  of  motion. 

Newmark's  method  was  used  to  integrate  the  finite 
element  equations  with  respect  to  time.  The  diagonal  mass 
matrix  was  used.  Two  types  of  impact  were  tried,  an 
inelastic  impact,  where  the  projectile  hits  the  plate  and 
sticks,  and  a Hertzian  impact,  where  the  force  between  the 
projectile  and  the  plate  is  governed  by  the  Hertzian  impact 
law.  Both  types  of  impact  produced  deflections  at  the 
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center  of  the  plate  which  agreed  well  with  deflections 
determined  experimentally  while  the  projectile  was  in 
contact  with  the  plate.  The  Hertzian  impact  method  allowed 
the  projectile  to  lose  contact  with  the  plate  and  impact  a 
second  and  possibly  a third  time.  A second  impact  has  been 
observed  experimentally. 

The  finite  element  program  was  verified  by  solving  the 
following  three  problems  which  have  analytical  solutions: 
static  deflection  of  an  isotropic  plate  with  a point  load 
at  the  center,  cylindrical  bending  of  an  isotropic  plate, 
and  free  vibration  of  an  isotropic  plate. 

The  finite  element  program  was  used  to  solve  plate 
impact  problems  which  do  not  have  analytical  solutions. 

The  deflections,  normal  stress,  and  average  shear  stress 
alone  the  .axis  resulting  from  impact  on  the  plate  were 
shown  graphically  for  isotropic,  3-layer,  and  15-layer 
composite  plates.  The  bending  stiffness  of  the  plate  and 
the  velocity  of  the  projectile  had  a large  effect  on  the 
disolacement  and  stresses.  The  effect  of  the  boundary 
conditions  on  the  displacements  and  stresses  increased  with 
time . 

The  normal  stresses  and  the  shear  stresses  show 
evidence  of  a large  flexural  wave  which  is  produced  by 
impact  of  the  plate  and  proDagates  outward.  The  flexural 
wave  may  be  the  cause  of  e xoerimenta 1 ly  observed 
delamination  crack  propagation. 
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CHAPTER  ONE 
INTRODUCTION 

Laminated  Composite  Plate 

A laminated  composite  material  consists  of  layers  of 
at  least  two  different  materials,  which  are  bonded 
together.  Lamination  is  used  to  combine  the  best  aspects 
of  each  layer  to  produce  a more  useful  material.  The 
properties  that  can  be  maximized  are  strength,  stiffness, 
low  weight,  impact  resistance,  etc.  The  layers  may  be 
anisotropic,  which  means  that  the  material  properties  are 
dependent  on  direction.  The  layers  can  be  chosen  to  tailor 
the  stiffness  and  strength  to  a specific  design  requirement 
of  a structure. 

A fiber-reinforced  composite  material  (FRCM)  consists 
of  fibers  in  a matrix.  If  the  fibers  are  laid  in  one 
direction,  then  the  material  is  anisotropic.  A laminated 
fiber-reinforced  composite  material  (FRCM)  plate  consists 
of  layers  of  FRCM  with  the  fiber  direction  in  different 
layers  oriented  in  different  directions,  which  gives  the 
layers  different  strengths  and  stiffness  in  the  various 
directions.  Fiber-reinforced  composites  can  be  designed  to 
have  the  advantage  of  high  strength-to-weight  ratio  and 
high  stif fness-to-weight  ratio,  and  can  also  be  designed  so 
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that  the  laminated  plate  may  have  preferential  directions 
of  enhanced  stiffness  and  strength.  Because  of  these 
advantages  FRCM  is  replacing  traditional  materials  in  many 
applications  in  aircraft,  automobiles,  and  sports 
equipment . 


Impact  on  a Laminated  Plate 

In  some  applications,  laminated  plates  made  of  fiber- 
reinforced  composite  material  (FRCM)  will  be  subject  to 
impact.  For  example,  the  leading  edge  of  an  aircraft  wing 
or  a jet  engine  blade  may  be  hit  by  foreign  objects  such  as 
stones,  nuts,  bolts,  hailstones  or  birds.  In  these  cases 
the  impact  resistance  of  layered  FRCM  plates  must  be  known 
in  order  to  develop  FRCM  plates  which  meet  aircraft  design 
safety  requirements. 

When  a moderately  thin  or  thick  plate  is  impacted  by  a 
projectile,  flexural  waves  are  produced.  Flexural  waves 
produce  displacements  and  stresses  in  the  plate  that  change 
with  time  and  are  different  from  the  stresses  produced  by  a 
static  load. 

In  this  study  a dynamic  finite  element  method  was  used 
to  calculate  the  displacements  and  stresses  in  a composite 
plate  impacted  by  a projectile.  It  was  assumed  that  the 
plate  remains  elastic,  and  that  no  failure  occurs.  Small 
deflection  plate  theory  is  used,  and  damping  in  the  plate 
is  assumed  to  be  negligible.  The  finite  element  stiffness 
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matrix  was  derived  from  the  laminated  anisotropic  plate 
theory  of  Yang,  Norris,  and  Stavsky  (1966),  which  includes 
shear  deformation  and  rotary  inertia.  Both  elastic  and 
inelastic  impacts  were  considered.  Newmark's  method  was 
used  to  integrate  with  respect  to  time. 

The  program  developed  has  the  capacity  to  take  into 
account  variations  in  the  material  properties  of  the 
laminas,  the  orientation  of  the  laminas,  the  number  of 
laminas,  the  density  of  the  material,  the  plate  dimensions, 
the  plate  boundary  conditions,  the  mass  of  the  projectile, 
the  velocity  of  the  projectile,  and  the  response  of  the 
material  to  indentation.  The  purpose  of  the  research 
project  was  to  calculate  the  stresses  in  a plate  impacted 
by  a projectile.  Shear  stresses  averaged  through  the 
thickness  and  normal  stresses  on  the  surface  were 
calculated  as  functions  of  time  and  position.  Knowledge  of 
where  and  when  the  stresses  are  a maximum  should  be  useful 
in  determining  where  failure  occurs  and  how  to  design 
composite  plates  with  maximum  resistance  to  impact. 

A Review  of  Work  on  Impact  of  Plates 

The  book  Impact  by  Goldsmith  (I960)  provides  a 
comprehensive  source  of  information  on  analytical  methods 
applied  to  impact  of  isotropic  plates.  A survey  of  work  on 
impact  of  composite  plates  is  given  by  Takeda  and 
Sierakowski  (1980). 
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Plate  impact  experiments  have  been  conducted  at  the 
University  of  Florida.  A small  steel  cylinder  is  fired 
from  an  air  gun  at  velocities  of  10-80  m/sec.  The  target 
is  a FROM  plate  as  shown  in  Figure  1.1.  The  plate  is 
clamped  at  the  edges.  The  deflection  and  surface  strains 
were  recorded  at  selected  points  in  some  of  the 
experiments.  At  subperforation  speeds  the  plate  usually 
fails  by  delamination  of  the  layers,  which  is  caused  by 
interlaminar  shear  and  normal  stresses.  In  the  flexural 
mode  of  deformation,  it  is  believed  that  the  shear  stress 
is  the  most  significant  cause  of  delamination. 

The  following  is  a survey  of  experimental  work  on 
impact  of  layered  composite  plates  conducted  at  the 
University  of  Florida.  The  failure  mechanisms  of  impacted 
layered  composite  plates  are  reported  by  Cristescu, 
Malvern,  and  Sierakowski  (1975)»  Sierakowski,  Malvern,  and 
Ross  (1976),  Ross,  Cristescu,  and  Sierakowski  (1976),  and 
Malvern,  Sierakowski,  and  Ross  (1979).  The  effect  of 
impactor  velocity  and  mass,  and  ply  orientation  on  plate 
delamination  is  reported  by  Takeda,  Sierakowski,  and 
Malvern  (1981a).  The  inplane  and  flexural  waves  produced 
by  impact  of  a composite  plate  are  reported  by  Takeda, 
Sierakowski,  and  Malvern  (1981b).  The  velocity  of  the 
delamination  crack  and  the  generator  strip  formation 
produced  by  impact  of  a composite  plate  is  reported  by 
Takeda,  Sierakowski,  Ross,  and  Malvern  (1982). 
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Figure  1.1.  Plate  Impacted  by  a Projectile 
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The  following  is  a survey  of  some  of  the  experimental 
work  on  impact  of  layered  composite  plate  excluding  work  at 
the  University  of  Florida.  Daniel,  Liber,  and  LaBedz 
(1979)  determined  experimentally  the  wave  velocities  in  a 
composite  plate  impacted  by  a high  speed,  low  density 
projectile.  They  determined  five  wave  velocities,  two  in- 
plane and  three  related  to  flexural  plate  deformation,  and 
compared  these  wave  velocities  with  theoretical  values. 

Sun,  Sankar,  and  Tan  (1981)  determined  the  contact  law 
between  a composite  plate  and  a steel  ball  for  static 
loading. 

The  dynamic  response  of  a layered  composite  plate  to 

<• 

impact  was  investigated  by  analytical  and  numerical 
methods.  Chow  (1971)  used  Laplace  transforms  to  analyze  a 
simply  supported  laminated  plate  subjected  to  a 
concentrated  impulse  load.  Moon  (1972)  studied  the  effect 
of  the  orientation  of  plies  in  each  layer  in  a composite 
material  on  the  extensional  and  bending  wave  velocities  and 
wave  fronts  produced  by  impact.  Kim  and  Moon  (1979)  used 
Laplace  transforms  in  time  and  Fourier  transforms  in  space 
to  analyze  stress  waves  in  an  anisotropic  plate  subjected 
to  a force  distributed  along  a line. 

Goldsmith  (I960)  used  the  normal  modes  method  to 
determine  the  dynamic  response  of  an  isotropic  plate  or 
beam  to  a rigid  impactor.  McQuillen,  Gause,  and  Llorens 
(1976)  extended  the  method  of  Goldsmith  (I960)  to  a 
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composite  beam  and  assumed  an  initial  displacement  and 
velocity.  Timoshenko  (1913)  used  normal  modes  and  a Hertz 
contact  law  to  analyze  the  deflections  of  a beam  impacted 
by  a projectile.  The  resulting  nonlinear  integral 
equations  were  solved  by  numerical  integration.  Sun  and 
Chattopadhyay  (1975)  extended  Timoshenko's  method  to  a 
laminated  simply  supported  composite  plate.  Dobyns  and 
Porter  (1981)  used  the  method  of  Sun  and  Chattopadhyay 
(1975)  for  a plate  subjected  to  a load  over  a small  area. 
Ramkumar  and  Chen  (1983)  used  Fourier  integral  transforms 
to  find  the  response  of  an  infinite  anisotropic  laminated 
plate  to  an  experimentally  determined  impact  force.  Keer 
and  Woo  (1984)  used  integral  equations  to  calculate  the 
pressure  distribution  in  an  isotropic  circular  plate 
impacted  by  a projectile  with  a large  radius  of  curvature 
producing  a large  impact  area.  Moyer  and  Gashaghai-Abdi 
(1984)  used  the  finite  difference  method  to  find  the 
dynamic  response  of  an  isotropic  plate  subject  to  a 
prescribed  initial  velocity.  Lee,  Du,  and  Liebowitz  (1984) 
analyzed  a composite  plate  subjected  to  an  assumed  pressure 
loading. 


Composition  of  the  Dissertation 
The  plate  equations  of  motion  of  a layered  composite 
material  are  derived  from  three  dimensional  elasticity 
theory  in  Chapter  Two.  Wave  motion  in  an  extended  medium 
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and  in  beams  and  plates  is  discussed  in  Chapter  Three.  The 
plate  equations  of  motion  are  derived  by  using  Hamilton's 
variational  principle  in  Chapter  Four,  and  extended  to 
finite  element  plate  equations  in  Chapter  Five.  A finite 
element  plate  equation  is  derived  by  Galerkin's  method  in 
Chapter  Six.  Newmark’s  method  of  integration  with  respect 
to  time  is  discussed  in  Chapter  Seven,  and  also  impact 
conditions.  The  steps  in  a finite  element  analysis  are 
given  in  Chapter  Eight.  Results  of  the  finite  element 
program  are  compared  with  known  solutions  in  Chapter 
Nine.  Displacements  and  stresses  are  presented  graphically 
in  Chapter  Ten  for  several  cases  of  impact  by  a projectile 
on  a composite  plate  and  compared  with  some  experimental 
results  previously  obtained  at  the  University  of  Florida. 


CHAPTER  TWO 

PLATE  EQUATIONS  OF  MOTION 

Shear  Deformable  Plate  Theory  With  Rotary  Inertia 

The  simplest  plate  theory  is  classical  plate  theory 
(CPT),  which  is  based  on  the  Kirchhoff  hypothesis  that 
normals  to  the  midsurface  before  deformation  remain 
straight  and  normal  after  deformation.  Classical  plate 
theory  underestimates  the  deflection  and  overpredicts  the 
natural  frequencies  because  transverse  shear  deformation 
and  rotary  inertia  are  neglected. 

The  omission  of  rotary  inertia  in  CPT  gives  the 
physically  impossible  result  that  a suddenly  applied  load 
produces  immediate  deflections  far  from  the  point  of 
application.  When  rotary  inertia  is  neglected,  flexural 
waves  of  infinitely  short  wave  length  propagate  at  an 
infinite  velocity,  which  is  physically  impossible. 

A shear  deformable  plate  theory  is  needed  when  the 
shear  deformation  is  significant.  Shear  deformation 
increases  as  the  ratio  of  thickness  to  length  of  the  plate 
increases.  Shear  deformation  is  larger  in  fiber  reinforced 
composite  material  (FRCM)  plates  than  in  isotropic  plates 
because  in  FRCM  the  ratio  of  Young's  modulus  to  shear 
modulus  can  be  very  large.  Shear  deformation  is  large  when 
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the  area  of  loading  is  small.  Shear  deformation  and  rotary 
inertia  are  important  for  impact  loading  because  waves  of 
short  wave  length  are  produced. 

Shear  deformable  plate  theory  is  based  on  the 
following  assumptions: 

(i)  the  deflections  are  small 

(ii)  normals  to  the  raidplane  before  deformation 
remain  straight  but  not  necessarily  normal  to 
the  midsurface  after  deformation 

(iii)  strains  normal  to  the  midplane  are  negligible. 
Assumption  (ii)  allows  the  plate  to  deform  in  transverse 
shear  as  well  as  bending  and  to  produce  shear  strains  as 
well  as  longitudinal  strains.  An  average  shear  stress 
through  the  thickness  can  be  calculated  from  the  shear 
strain . 

Books  by  the  following  authors  are  useful  sources  of 
information  on  isotropic  and  anisotropic  plates: 

Timoshenko  and  Woinowsky-Krieger  (1959),  Leissa  (1969),  and 
Washizu  (1975).  Books  and  articles  by  the  following 
authors  are  useful  sources  of  information  on  layered 
anisotropic  plates:  Ashton  and  Whitney  (1970),  Jones 

(1974),  Vinson  and  Chou  (1975),  and  Ross,  Sierakowski,  and 
Sun  (1980).  A review  of  the  development  of  plate  equations 
for  layered  anisotropic  materials  is  presented  by  Vinson 
and  Chou  (1975,  p.  228)  and  Reddy  (1980). 
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Timoshenko  (1921)  beam  theory,  which  includes 
transverse  shear  and  rotary  inertia  effects,  was  extended 
to  isotropic  homogeneous  plates  by  Reissner  (1945)  and 
Mindlin  (1951).  This  shear  deformable  plate  theory  was 
extended  to  laminated  anisotropic  plates  by  Yang,  Norris, 
and  Stavsky  (1966)  and  is  referred  to  as  YNS  theory. 
Ebcioglu  (1964)  extended  Von  Karman  plate  theory  for 
isotropic  plates  to  include  the  effects  of  rotary  inertia 
and  shear  deformation  in  the  theory  of  anisotropic  plates. 

The  partial  differential  equations  based  on  YNS  theory 
can  be  solved  in  closed  form  for  only  a few  special 
cases.  Whitney  and  Pagano  (1970)  have  solved  plate 
equations  for  a static  sinusoidally  distributed  load  and 
for  free  vibration  of  a simply  supported  rectangular  plate. 

The  dynamic  response  of  an  infinitely  long  simply 
supported,  laminated  composite  plate  in  cylindrical  bending 
was  analyzed  under  the  following  loading  conditions:  Sun 

and  Whitney  (1974),  and  Sun,  Whitney,  and  Whitford  (1975) 
used  a suddenly  applied  dynamic  pressure.  Sun  and  Whitney 
(1976)  used  an  initial  inplane  stress  and  a suddenly 
applied  concentrated  load.  Whitney  and  Sun  (1977)  used  a 
uniform  or  line  load  applied  as  a pulse. 

Finite-difference  methods  have  been  used  to  solve  the 
equations  of  classical  plate  theory  (CPT)  for  anisotropic 
laminated  plates.  Sierakowski,  Chang-Tsan  Sun,  and  P.W. 

Sun  (1973)  computed  the  maximum  deflection  of  a layered 
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plate  under  a static  concentrated  load.  Johnson  and  Woolf 
(1984)  computed  the  static  deflections  and  stresses  in  a 
thin  orthotropic  plate  under  uniform  and  point  loads. 

Much  work  has  been  done  to  develop  plate  finite 
elements  based  on  classical  plate  theory  (CPT),  shear 
deformable  plate  theory,  and  Yang,  Norris,  and  Stavsky 
theory  (YNS).  Shear  deformable  plate  theory  and  YNS  theory 
are  better  than  classical  plate  theory  (CPT)  for  finite 
elements  because  the  continuity  requirement  is  reduced  from 
to  Cq.  A continuity  requirement  of  Cq  requires  that 
only  the  displacements  and  not  the  derivatives  be 
continuous.  A review  of  plate  bending  elements  is  given  by 
Hrabok  and  Hrudey  (1984). 

A review  of  work  done  to  develop  plate  finite  elements 
based  on  Yang,  Norris,  and  Stavsky  theory  (YNS)  is  provided 
by  Lakshminarayana  and  Murthy  (1984).  Pryor  and  Barker 
(1971)  used  a finite  element  with  7 degrees  of  freedom  at 
each  node.  Reddy  (1980)  presented  a penalty  plate  bending 
element  based  on  YNS  theory. 

Plate  Strain-Displacement  Relationships 

Consider  a plate  of  constant  thickness,  h,  whose  edge 
view  is  shown  in  Figure  2.1.  The  x and  y coordinate  axes 
lie  in  the  middle  plane  of  the  undeformed  flat  plate.  In 
the  undeformed  position,  the  z axis  is  perpendicular  to  the 
plate  as  shown  in  Figure  2.1.  In  the  deformed  position,  a 
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Edge  View  of  a Plate  in  the  Initial  Unstressed 
Position  and  in  the  Final  Position  with  Shear 
and  Bending  Deformation 


Figue  2.1. 
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point  initially  on  the  midplane  has  moved  a distance  u°  in 
the  x direction  and  w in  the  z direction.  The  displacement 
v°  in  the  y direction  is  not  shown.  In  the  deformed 
position,  the  line  initially  perpendicular  to  the  x axis 
has  rotated  an  angle  4>x  measured  from  the  z axis  toward  the 
x axis. 

A point  A which  is  located  a distance  z above  the  mid- 
surface  in  the  undeformed  plate  moves  to  point  A'  in  the 
deformed  plate.  The  x-component  displacement  of  point  A is 

u ( x, y , z , t ) = u° ( x , y , t ) + z sin  *x(x,y,t)  2.1 

For  small  deformations  we  assume 

sin  ij;  = i}»  2.2 

x x 


Then 


u(x,y,z,t)  = u° ( x , y , t ) + z 4’x(x,y,t)  2.3 

Similarly  by  considering  a point  in  the  yz  plane  we  obtain 
v ( x , y , z , t ) = v° ( x , y , t ) + z <|>  (x,y,t)  2.4 

The  transverse  deflection  w is  assumed  to  be  constant 


through  the  thickness. 
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v(x,y,z,t)  = w°(x,y,t)  = w(x,y,t)  2.5 

where  u°,  v°,  and  w°  are  the  midplane  displacement 
components  in  the  xf  y,  and  z directions,  respectively,  and 
<px  and  4»y  are  the  rotations  of  the  normals  in  the  x and  y 
directions,  respectively. 

The  rotation  ij;x  is  the  sum  of  the  slope  of  the  plate, 
-3w/3x,  and  the  shear  deformation  yxz*  The  negative  sign 
is  included  because  the  slope  3w/3x  is  negative  when  the 
rotation  is  positive. 


'i' 


x 


3w 

3x  + 


xz 


2.6a 


3 w 

ill  = - — — + y 

y 9y  yz 


2.6b 


In  Equation  2.6  we  have  assumed  that  u and  v vary  linearly 
through  the  thickness  direction  while  w is  constant  through 
the  thickness. 

The  displacement  Equations  2.3,  2.4,  and  2.5  are 
substituted  into  the  definition  of  linear  strain  to  obtain 
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yz 


aw 

ay 


Since  w,  $x,  and  ipy  are  independent  of  z,  the  engineering 
shear  strains  yxz  and  YyZ  are  constant  through  the 
thickness.  The  assumption  that  w is  independent  of  z 
implies  that  ez  is  zero.  The  equations  presented  above  are 
based  on  the  original  development  by  Reissner  (1945). 


For  the  constitutive  equation  of  a lamina  and  a 
laminate  presented  in  this  and  the  following  subsection, 
see,  for  example,  Ross,  Sierakowski,  and  Sun  (1980), 

Chapter  Three  by  Sun.  Consider  a unidirectionally 
reinforced  lamina  with  the  axes  x^ , x?,  x ^ located  in  the 
principal  material  directions  of  the  orthotropic  lamina  as 
shown  in  Figure  2.2.  The  principal  material  directions  are 
parallel  to  the  intersections  of  the  three  mutually 
orthogonal  planes  of  material  symmetry.  For  a 
unidirectionally  fiber- reinforced  composite  material,  the 
three  principal  material  directions  are  chosen  as 
follows:  x-]  in  the  direction  of  the  fiber,  X2 

perpendicular  to  the  fiber  and  parallel  to  the  plate,  and 
x^  perpendicular  to  the  fiber  and  the  plate.  The  x>j-X2 
plane  is  chosen  to  coincide  with  the  midplane  of  the 


Constitutive  Equations  for  a Lamina 


lamina. 
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The  constitutive  equation  for  a lamina  is 


The  normal  strain  is  zero,  and  a ^ is  omitted.  The 
stresses  are  shown  in  Figure  2.3.  The  reduced  stiffnesses, 
Qj^,  are  the  reduced  elastic  constants  in  the  principal 
material  directions. 


Q11 

“ E1/(1  - v12v21) 

2.9a 

Q22 

= E2/(1  - »12v21) 

2.9b 

Q12 

= v-12E2'/('  1 ' v12v2l'  = v21E1^1  " v12v21* 

2.9c 

Q44 

' G23 

2.9d 

Q55 

' °13 

2.9e 

Q66 

- G12 

2.9f 

where  and  E^  are  Young's  moduli  in  the  and  x? 

directions,  respectively.  Poisson's  ratio  v.  . is  for 

* d 

transverse  strain  in  the  x . direction  when  stressed  in  the 

J 
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Figure  2.3.  Stresses  in  a Lamina 
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x^  direction,  and  G^,  3nd  are  the  shear  moduli  in 

the  X2~x^,  x-j-Xtj,  and  x-j-X2  planes,  respectively. 

The  principal  material  axes  of  the  lamina  (x-j,X2»x^) 
do  not  in  general  correspond  to  the  reference  axes 
(x,y,z).  The  stress-strain  relation  for  an  arbitrarily 
oriented  lamina,  shown  in  Figure  2.4,  is  given  by 


where  are  the  transformed  reduced  stiffnesses.  The 

stresses  are  shown  in  Figure  2.9. 

The  transformed  reduced  stiffnesses  Q—  with  respect 
to  a coordinate  system  x,  y,  and  z,  rotated  an  angle  9 from 
the  principal  material  axes  are  as  follows: 

0^  = Q1 1cos49+2(Q12+2Q-^)sin20cos29+Q22sin40 
Q1?  = (Q11+Q2p-4Q6g)sin20cos20+O1?(sin40+cos40) 

Q?2  = Q11sin4e+2(Q1p+2Qgg)sin29cos20+Q?9cos'40 
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§16  = (Q-11  -Q12'2<366^sinecos^e+^Q12-Q22+2Q66^sin39cos0 

§26  = (^11 ~Qi2~-^66^  s^n^0cos9  + ^i2-^22+^^66^  s^n0cos^9  2,11 
§44  = Q44cos2e+Q55sin2e 

§45  = (Q44-Q55)cosesine 

§55  = Q55cos29+Q44sin2e 

0 0 A A 

§66  = (Qii+Q22"2!^12"2<^66^sin  0COS  0+(^66^sin  9+cos40^ 

Constitutive  Equations  for  a Laminate 
In  the  previous  section  stress-strain  relationships 
for  a single  lamina  were  discussed.  A laminate  consists  of 
two  or  more  laminas  bonded  together  to  form  a plate.  The 
principal  material  axes  of  different  laminas  may  be 
oriented  in  different  directions.  The  properties  of  the 
laminate  depend  on  the  properties  of  each  lamina  and  on  the 
orientations  of  the  various  laminas  in  the  plate.  In  this 
section  constitutive  equations  for  a laminate  will  be 
discussed . 

As  an  example,  consider  the  regular  symmetric  cross- 
ply  laminate  shown  in  Figure  2.6.  It  is  composed  of  three 
identical  orthotropic  laminas  with  the  center  lamina 
rotated  an  angle  of  90°.  It  is  symmetric  with  respect  to 
the  midplane.  For  the  top  and  bottom  laminas,  where 
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Figure  2.4.  Arbitrarily  Orientated  Lamina 


z 


y 


Figure  2.5.  Stresses  in  a Laminate 
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*2 


Figure  2.6.  Three-Layer  Cross-ply  Symmetric  Laminate 
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rotation  is  0°,  the  transformed  reduced  stiffness  is 

the  same  as  the  reduced  stiffness  For  the  center 

lamina  with  a rotation  of  90°, 


2.12 


The  stress-strain  relationship  for  a single  lamina, 
Equation  2.10,  remains  valid  for  a lamina  in  a laminate, 
and  combined  with  the  strain  displacement  relationships. 
Equation  2.7,  yields  for  the  nth  layer 


The  stresses  in  the  nth  layer  are  defined  in  terms  of  the 
in-plane  displacements,  u°,  v°,  the  out-of-plane 
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displacement,  w,  and  the  rotations,  <j>x  and  ij>y.  Note  that 
these  quantities  are  assumed  constant  throughout  the 
thickness  of  the  laminate.  Therefore  the  strains  are 
linear  in  z,  but  the  stresses  may  be  discontinous  between 
laminas  because  the  Q^^n^may  be  different  for  different 
laminas . 

The  force  and  moment  resultants  acting  on  the  laminate 

are  found  by  integrating  the  stresses,  given  for  all 

laminas  by  Equation  2.13,  through  the  thickness  of  the 

laminate,  h.  The  tractions  Nv,  N„  and  Nv„  with  dimension 

x y xy 

(F/L)  are  defined  as 


The  transverse  shearing  forces  Qx  and  Qy  with  dimension 
(F/L)  are  defined  as 


The  bending  moments  Mx  and  My  and  twisting  moment  MXy  with 
dimension  (FL/L  = F)  are  defined  as 
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2.16 


The  resultants  are  shown  in  Figure  2.7. 

The  stress-strain  relation  for  the  nth  lamina, 
Equation  2.13,  is  substituted  into  the  definitions  of  the 
resultants,  Equations  2.14,  2.15  and  2.16,  to  obtain 
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The  plate  stiffnesses  A.  B.  . and  D.  . are  defined  as 

ij  10  ij 
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Figure  2.7.  Force  and  Moment  Resultants  Acting  on  a Plate 

a)  Normal  and  Shear  Resultants 

b)  Moment  Resultants 
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(i,0=4,5)  2.18b 


constants  of  the  nth  lamina  in  the  plate  coordinates  and  < 
is  the  shear  correction  factor.  Reissner  (1945) 


plates.  Whitney  and  Pagano  (1970)  used  a value  of  5/6  for 
layered  composite  plates.  A value  of  5/6  was  used  in  this 
study.  The  bounds  of  integration  zn  and  zn_ ^ are  defined 
as  the  z-coordinates  measured  from  the  midplane  of  the 
laminate  to  the  surfaces  of  the  nth  lamina  as  shown  in 
Figure  2.8.  After  integration.  Equations  2.18  become 


recommended  a value  for  < ^ of  5/6  for  homogeneous  isotropic 


(i, ,3  = 1 ,2,6) 


(i, 0=1, 2, 6) 


2.19 


(i, 0 = 1, 2, 6) 
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(i, 0=4,5) 


layer  number 


Figure  2.8.  Edge  View  of  an  N Layer  Laminate 
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The  are  called  the  extensional  stiffnesses,  the  are 

called  the  coupling  stiffnesses,  and  the  are  called  the 


bending 

stiffnesses . 

Equations  of  Motion 

The  two-dimensional  equations  of  motion  of  a plate  are 


derived 

from  the  three-dimensional  equations  of  motion: 

o + t + t = pu  , , 2.20a 

x,x  xy,y  xz,z  ,tt 

t + a + t = pv  . , 2.20b 

xy, x  y , y yz,z  ,tt 

t + t + a = pw  , , 2.20c 

xz, x  yz,y  z,z  ,tt 

A* 

where  p 

is  the  density,  mass  per  unit  volume.  The  comma  is 

used  to 

denote  differentiation  with  respect  to  x,  y,  z or 

t. 

Three  of  the  five  plate  equations  of  motion  are 
derived  by  integrating  Equations  2.20  with  respect  to  z and 
using  the  resultant  definitions,  Equations  2.14  and  2.15. 


The  two 

inplane  equations  of  motion  are 
Nx,x  * Nxy,y  ■ mu,tt  * R*x,tt  2-21 

Nxy,x  * Ny,y  * mv,tt  + My,tt  2'22 

2.22 
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where  mass  per  unit  area  is 


h/2 

m = J P dz 

-h/2 


For  uniform  density 


m = Ph 


One  out-of-plane  equation  of  motion  is 


Q + Q + q = mw 
x,x  y,y  ,tt 


2.23 


2.24 


2.25 


where  q2  is  the  applied  load,  a force  per  unit  area. 

The  other  two  out-of-plane  equations  are  derived  by 
multiplying  Equations  2.20  a and  b by  z,  integrating  with 
respect  to  z,  and  using  the  stress  resultant  definitions. 


2.26 

2.27 

The  rotary  interia  per  unit  area  is 


M +M  -Q  =1’!'  , , + Ru 

x,x  xy,y  x x,tt  ,tt 


M + M - Q = It  ,,  + Rv  , 

xy,x  y,y  *y  y,tt  ,tt 


h/2 

I = I 

-h/2 


2 

P z dz 


2.28 
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and  for  uniform  density 


2.29 


The  rotational  coupling  inertia  is 


h/2 


R = / 


p z dz 


2.30 


-h/2 


If  the  laminate  is  symmetric  with  respect  to  geometry 
and  the  material  properties  about  the  middle  surface  and 
the  forces  Nx,  NXy,  and  Ny  are  linear  functions  of 
displacement,  then  the  two  inplane  Equations  2.21  and  2.22 
are  uncoupled  from  the  three  out-of-plane  Equations  2.24, 
2.25  and  2.26.  The  reason  is  as  follows.  For  a symmetric 
laminate  the  coupling  stiffnesses,  B^,  are  zero.  Then 
according  to  Equation  2.1?  the  normal  resultants  are 
functions  of  u°  and  v°  only  and  not  of  w,  i|>x  or  <J>y.  This 
uncouples  the  equations. 

Classical  plate  theory  (CPT)  is  a special  case  of 
shear  deformable  plate  theory.  For  CPT  the  shear  strains, 
Yxz  and  Yyz,  are  zero  and  Equations  2.6  become 


2.31 
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and  the  displacements  are 

u(x ,y , 2, t ) = u° ( x , y , t ) 
v(x,y,z,t)  = v°(x,y,t) 
w(x,y,z,t)  = w°(x,y,t) 


- z -57  (x,y,t)  2.32a 

3 w 

~ z ly  2.32b 

2.32c 


CHAPTER  THREE 
WAVE  MOTION  IN  SOLIDS 

Types  of  Waves 

Transverse  impact  of  a plate  by  a projectile  produces 
transverse  vibrations  of  the  plate.  The  transverse 
vibrations  may  be  considered  as  a sum  of  mode  shapes  or  as 
flexural  wave  motion.  The  mode  shapes  are  easiest  to 
understand  but  are  most  useful  for  explaining  cyclic 
vibration.  Flexural  waves  are  the  best  way  to  explain 
transient  motion  of  a plate,  but  flexural  waves  are 
difficult  to  analyze  because  they  do  not  travel  with 
constant  shape  or  velocity.  In  order  to  understand  wave 
motion  in  a plate,  the  simpler  waves  in  an  extended  medium 
and  in  a beam  will  be  considered  first. 

Wave  Motion  in  an  Extended  Medium 

In  an  extended  isotropic  elastic  medium,  the  three- 

dimensional  equations  of  motion,  Equations  2.20,  have  a 

solution  in  terms  of  two  types  of  waves,  a dilatational  and 

an  equivoluminal  wave  according  to  Malvern  (1969,  p.  550). 

The  dilatational  wave  involves  no  rotation  and  travels  with 

1b 

velocity  [(i  + 2G)/p]  ^ . A dilatational  wave  in  air  is 

called  a sound  wave.  The  equivoluminal  wave,  also  called 
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a shear  wave,  distortional  wave,  or  rotational  wave, 
travels  at  velocity  (G/p)  ^2  . in  an  unbounded  isotropic 
elastic  solid  these  two  waves  are  the  only  two  possible 
waves, and  any  plane  wave  must  propagate  at  one  or  the  other 
of  these  two  velocities. 

Bounding  free  surfaces  affect  the  wave  motion  in  an 
elastic  medium.  If  there  is  one  plane  bounding  free 
surface,  then  there  is  also  a surface  wave  called  a 
Rayleigh  wave.  If  the  elastic  medium  is  a plate  with  two 
parallel  bounding  free  surfaces  closely  spaced  compared  to 
the  other  dimensions,  then  the  flexural  wave,  also  called 
the  bending  wave,  produces  larger  stresses  than  the 
dilatational  or  shear  wave.  Kim  and  Moon  used  analytical  " 
methods  to  study  waves  in  a layered  plate  and  concluded, 

"A  relatively  insignificant  part  of  the  stress  propagates 
via  dilatational  and  shear  waves.  The  maximum  stress  is 
transmitted  by  a bending  wave  which  propagates  in  the  plate 
at  a much  slower  speed"  (1979.  p.  1130). 

Wave  Motion  in  a Beam 

Wave  motion  in  a beam  is  similar  to  wave  motion  in  a 
plate  because  a beam  is  a plate  with  a large  length  to 
width  ratio.  Beam  theory  is  simpler  than  plate  theory 
because  a beam  has  only  one  dimension,  length,  whereas  a 
plate  has  two  dimensions,  length  and  width.  Wave  motion  in 
a beam  is  easier  to  understand  than  wave  motion  in  a plate 
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because  there  are  some  closed  form  solutions  for  particular 
geometries  and  loadings.  An  understanding  of  waves  in  a 
beam  will  be  useful  in  understanding  the  wave  motion  in  a 
plate  presented  in  Chapter  Ten. 

There  are  closed  form  solutions  for  the  equations 
governing  torsional  and  longitudinal  waves  in  beams,  and 
these  solutions  are  discussed  in  many  books  on  vibration. 
There  are  very  few  closed-form  solutions  for  equations 
governing  transverse  waves  in  beams.  Longitudinal  and 
flexural  waves  in  a beam  will  be  discussed  here  because 
these  same  types  of  waves  also  propagate  in  a plate. 

A longitudinal  wave  in  a beam  is  governed  by  the 
equation 


3.1 


The  solution  of  Equation  3.1  for  an  infinite  beam  is  given 
by  Berg  and  McGregor  (1966,  p.  206-209).  For  an  infinite 
beam  initially  at  rest,  with  initial  conditions 

u(x,0)  = f(x)  u,t(x,0)  = 0 3.2 


— oo 


< x < 


oo 
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the  general  solution  is 

f(x+c  t)  + f(x-c  t) 

u(x,t)  = ^ 2 — 3.3 

If  f(x)  is  a function  that  vanishes  outside  a finite 
interval,  then  Equation  3.3  represents  two  longitudinal 
waves  traveling  in  opposite  directions,  each  with  constant 
shape  defined  by  the  functions  f and  constant  velocity 


c = /E/p 
o 


3.4 


According  to  elementary  beam  theory,  which  neglects 
rotary  inertia  and  shear  deformation,  flexural  waves  in  a 
thin  beam  are  governed  by  the  equation 


3 


El 

pA 


34i 
3 x^ 


3.5 


Flexural  waves  are  governed  by  a partial  differential 
equation  (PDE),  which  is  4th  order  in  space  (Equation  3*5), 
whereas  longitudinal  and  torsional  waves  are  governed  by 
partial  differential  equations  which  are  2nd  order  in  space 
(Equation  3.1).  A PDE  which  is  4th  order  in  x does  not 
have  a solution  of  a wave  traveling  with  constant  amplitude 
and  shape  because  a 4th  order  PDE  does  not  have  solution 
like  Equation  3.3  and  EI/pA  does  not  have  the  dimensions  of 
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velocity.  The  solution  to  the  PDE  which  is  4th  order  in 
space  is  a wave  which  changes  shape  and  velocity  as  it 
propagates. 

A closed  form  solution  of  Equation  3«5  exists  for 
sinusoidal  flexural  waves  propagating  along  a beam.  The 
wave  velocity  or  phase  velocity  is  inversely  proportional 
to  the  wavelength.  This  means  that  waves  with  short 
wavelengths  travel  faster  than  waves  with  longer 
wavelengths. 

A flexural  wave  is  dispersive.  A flexural  pulse  may 
be  represented  by  a Fourier  integral  sum  of  sinusoids  of 
all  wave  lengths  and  different  amplitudes,  each  propagating 
at  its  own  velocity.  As  the  pulse  propagates,  it  changes 
shape  because  the  different  Fourier  components  travel  at 
different  velocities. 

Elementary  beam  theory  predicts  that  infinitely  short 
wavelength  sinusoidal  waves  would  propagate  at  an  infinite 
velocity,  which  is  physically  impossible.  Timoshenko 
(1921)  beam  theory,  which  includes  rotary  inertia  and  shear 
deformation,  correctly  predicts  the  wave  velocities  for 
sinusoidal  flexural  waves,  even  for  waves  with  wave  lengths 
as  small  as  the  thickness  of  the  beam.  Very  short  waves 
propagate  at  the  Rayleigh  wave  velocity,  which  is  close  to 
the  shear  wave  velocity  of  an  unbounded  medium. 
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Wave  Motion  in  a Plate 

In  a plate  governed  by  plate  equations,  there  are  two 
elastic  waves,  a longitudinal  wave  and  a flexural  wave. 
These  waves  are  produced  by  the  interaction  of  the  two 
three-dimensional  waves  with  the  boundary.  The 
longitudinal  wave  in  a plate  is  similar  to  the  longitudinal 
wave  in  a bar  and  to  the  dilatational  wave  in  an  extended 
elastic  medium.  The  velocity  of  the  longitudinal  wave  in 
an  isotropic  elastic  infinite  plate  is 


c 


L 


E 

p(l-v2) 


3.6 


provided  the  wave  length  is  long  compared  to  the 
thickness.  If  the  wavelength  is  very  small  compared  to  the 
thickness  of  the  plate  the  longitudinal  wave  travels  at  the 
Rayleigh  wave  velocity.  For  wavelengths  comparable  to  the 
thickness  of  the  plate,  dispersion  occurs  and  the  velocity 
depends  on  the  ratio  of  the  wavelength  to  thickness 
according  to  Kolsky  (1963,  p.  83).  The  longitudinal  wave 
velocity  in  a plate  differs  from  the  longitudinal  wave 
velocity  in  a beam  by  the  factor  (1-v^)”  ^2  , as  can  be  seen 
by  comparing  Equations  3.6  and  3.4. 

Longitudinal  waves  are  governed  by  the  two  inplane 
Equations  2.21  and  2.22  in  the  x and  y directions 
respectively.  For  a symmetric  layup  laminate  the  inplane 
Equations  2.21  and  2.22  are  uncoupled  from  the  out-of-plane 
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Equations  2.25,  2.26,  and  2.27.  Longitudinal  waves  are  not 
significant  for  transverse  impact  loading  of  a glass  epoxy 
laminate  with  a length  of  side  to  thickness  ratio  of  42, 
according  to  Takeda,  Sierakowski,  and  Malvern  (1981b). 

They  used  strain  gages  to  record  the  strain  on  the  surface 
of  a symmetric  layup  plate  caused  by  transverse  impact  of  a 
plate.  They  were  able  to  identify  strain  produced  by 
longitudinal  and  flexural  waves  by  identifying  their 
velocities.  The  strain  produced  by  longitudinal  waves  was 
very  small  compared  to  the  strain  produced  by  flexural 
waves.  Daniel,  Liber,  and  LaBedz  (1979)  also  agree  with 
this  conclusion.  They  determined  experimentally  that  the 
amplitude  of  the  inplane  wave  is  less  than  10  percent  of 
the  amplitude  of  the  flexural  wave. 

The  flexural  wave  is  the  predominant  wave  produced  by 
transverse  impact  of  a moderately  thin  or  thick  plate.  The 
three  out-of-plane  Equations  2.25,  2.26,  and  2.27  govern 
flexural  waves.  If  shear  deformation  is  neglected,  the 
three  out-of-plane  equations  reduce  to  one  equation 


M + 2M  + M 

x , xx  xy , xy  y,yy 


mw 


tt 


3.7 


and  in  terms  of  displacement  for  an  isotropic  medium 


DV  w +mw 


, tt 


= q. 


3.8 
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where 


V w = w + 2w  + w 

» xxxx  , xxyy  ,yyyy 

and  D is  the  bending  stiffness  of  the  plate 


3.9 


D 


Eh' 


12( 1-v  ) 


3.10 


Equation  3.8  is  based  on  classical  plate  theory,  which  is 
valid  for  thin  plates.  Equation  3.5  which  governs  beams  is 
a special  case  of  the  plate  equations  3.8  and  3.9,  where  El 
replaces  D and  the  y derivatives  are  omitted.  The  results 
of  finite  element  analysis  presented  in  Chapter  Ten  will 
show  flexural  waves  in  layered  composite  plates. 


CHAPTER  FOUR 

A VARIATIONAL  FORMULATION 
OF  THE  PLATE  EQUATION  OF  MOTION 


Hamilton’s  Principle 

In  Chapter  Two  the  plate  equations  of  motion  were 
derived  from  the  three-dimensional  equations  of  motion.  In 
this  chapter  the  plate  equations  of  motion  will  be  derived 
from  Hamilton's  principle,  which  is  a variational 
principle.  Both  a variational  formulation  and  a 
differential  formulation  of  the  plate  equations  of  motion 
will  be  derived.  The  differential  formulation  of  the  plate 
equations  of  motion  will  be  the  same  as  the  plate  equations 
of  motion  in  Chapter  Two  and  will  serve  as  a check  on  the 
algebra.  The  variational  formulation  of  the  plate 
equations  of  motion  will  be  used  in  Chapter  Five  to  derive 
the  finite  element  stiffness  and  mass  matrices. 

Hamilton's  principle,  as  stated  by  Langhaar,  is  as 
follows : 

Among  all  motions  that  will  carry  a conservative 
system  from  a given  initial  configuration  Xq  to  a 
given  final  configuration  in  a given  time 
interval  (tQ,  t>|),  that  which  actually  occurs 
provides  a stationary  value  to  the  integral. 

(1962,  p.  239) 
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/ 1 (6T  - 6U  - 6W)dt  = 0 4.1 

*0 

where  T is  the  kinetic  energy  of  the  body,  U is  the  strain 
energy  of  the  body,  and  W is  the  work  done  by  the  body. 

The  variational  symbol,  5,  means  the  first  variation  of 
what  follows. 

Strain  Energy  in  a Plate 

The  total  internal  strain  energy  of  a linear  elastic 
body  is 


U 


4.2 


The  first  variation  of  the  strain  energy  for  a linear 
elastic  body  is 


<5(J  = 


4.3 


Since 


e 

z 


= 0 


4.4 


then 
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h/2 

6U  = / / / ( a fie  + ct  fie  + t fiy 

-h/2  x x y y xy  'xy 


+ t fiy  + t fiy  ) dxdydz 
xz  xz  yz  yz 


4.3 


where  h is  the  thickness  of  the  plate.  We  use  the  strain 
displacement  relationships  (Equation  2.7)  to  obtain 


h/2 

5U  = / ff(ax(6u°y  + Z'S'l'y  y ) + + Z 5 * ) 

hf  2 x >x  x»x  y >y  y»y 


+ T [ 6u°  + fiv°  + Z ( fi  it  + 5<J>  )] 

xy  , y ,x  x , y ry,x 


4 • 6 


+ tyv6^w  y + *Y)  + t 6(v  v + 'I'v ) jdxdydz 
aj  9 A a y ^ »y  y 


The  first  variation  of  the  strain  in  terms  of  the  stress 
and  moment  resultants  defined  in  Equations  2.14,  2.15  and 
2.16  is 


<5U  = // [N  <5u°  + N 5v°  + N 5(u°  + v°  ) 

X ,x  y ,y  xy  , y ,x' 


+ M5'4'  + M fiiji  +M  5 ( +i{)  ) 

x x,x  y y,y  xy  KYx,y  y,x 


4.7 


+ + w ) + Q fi(*  + w ) ]dxdy 

A A f A y y f y 


Equation  4.7  is  the  variational  form  of  the  plate  equations 
of  motion  and  will  be  used  in  Chapter  Five  to  derive  the 
finite  element  equations. 


The  two  integration-by-parts  formulas  for  double 
integrals,  as  stated  by  Zienkiewicz  (1977,  p.  94),  are 
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4.9 


4.8 


where  nx  and  ny  are  the  components  (direction  cosines)  of 
the  unit  normal 


on  the  boundary  and  ds  is  the  arc  length  of  an 
infinitesimal  segment  of  the  boundary.  The  integral  is 
taken  in  an  anticlockwise  direction  around  the  boundary. 
Equations  4.8  and  4.9  were  applied  to  Equation  4.7  to 
obtain 


n = n i + n j 
x y° 


4.10 


(4.11  ) 


+ ^[(Nxnx  + Nxyny)5u 


+ Q n ) <5w  +(Mn  + M n )<5<l> 
y y XX  xy  y;  X 


+ (M  n + M n )8<(i  Ids 
xy  x y y y 
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Work  Done  on  a Plate 

The  external  virtual  work  done  on  a plate  is  the  sum 
of  the  virtual  work  done  by  forces  acting  through  a virtual 
displacement  and  the  virtual  work  done  by  moments  acting 
through  a virtual  rotation,  with  the  forces  and  moments 
held  constant  during  the  infinitesimal  virtual  motion. 

Since  the  stress  resultants,  N,  Q,  and  M,are  defined  as 
forces  or  moments  per  unit  length  along  the  edge  of  the 
plate,  the  stress  resultants  are  integrated  over  the  length 
of  the  plate  on  which  the  work  is  done.  The  surface  load, 
qz , with  units  of  force  per  unit  area  is  integrated  over 
the  plate  area. 

The  virtual  work  done  on  the  body  by  external  loads, 
stress  resultants,  and  moment  resultants  is 

6W  = - <f[(N’n  + N'  n )6u°  + ( N T n + N'n  )6v° 

JLVxx  xy  y xy  x y y' 

+ (Q'n  + Q'n  )<5w  + (M’n  + M*  n )5<|/ 

^x  x y y xx  xy  y x 

+ (M'  n + M'n  )6<J»  ]ds  - //  q 6w  dxdy  4.12 

xy  x y y y J 1 z J 

where  the  primed  quantities  are  prescribed  edge  loads. 


46 


Kinetic  Energy  of  a Plate 
The  kinetic  energy,  T,  of  a body  is 

t = rrr  £ ii  . dV 

JJJ  2 at  at 


4.13 


where  p is  the  density  and  u is  the  displacement  vector 


AAA 

u = ui  + vj  + wk 


4.14 


Hamilton's  principle  requires  that  we  find  the  first 
variation  of  the  kinetic  energy  integrated  from  initial 
time  tg  to  a final  time  t-j 

6/  1 T dt  = / 1 [fffp  -§£  • 5(i^)dV]dt  4.15 

fc0  t0 

Integration  by  parts  yields 


5/  1 Tdt  = ///p[-|t-  * 6u 
t,  z 


' 0 


fc1  _ jt1 

tQ  ' tQ  3t2 


• 5u  dt]dV 


4.16 


There  is  no  variation  in  the  displacement  at  the  initial 
and  final  times  because  Hamilton's  principle  assumes  that 
the  initial  and  final  configurations  are  prescribed. 


6u(tQ)  = 5u( t1 ) = 0 


4.17 
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Therefore , 


t 

«J 


1 


Tdt 


///[p/ 


32a 


*0  at4 


Su  dt ]dV 


4.18 


Equation  4.18  can  be  expressed  in  terms  of  the 
displacements,  u,  v,  and  w,  of  Equation  4.14,  the  mass,  m, 
of  Equation  2.23  and  the  rotary  inertia,  I,  of  Equation 
2.28,  and  the  rotational  coupling  inertia,  R,  of  Equation 
2.30. 


t1  t1 

5/  Tdt  = - / //[ m(u  + v ff<5v°  + w t-t-<5w) 


, tt 


, tt 


* R(*x,ttSu°  + *y,tt5',°  * u,tt5*x  + v,tt6^y)  4-19 


* I(*x,ttS*x  * *y,ttS,,y)ldAdt 


Plate  Equations  of  Motion 
Hamilton's  principle  is  now  used  to  find  the 
variational  form  of  the  plate  equations  of  motion.  The 
first  variation  of  the  strain  energy,  Equation  4.11,  the 
kinetic  energy.  Equation  4.19  and  the  work,  Equation  4.12, 
are  substituted  into  Hamilton's  principle,  Equation  4.1, 
and  the  coefficients  of  the  variational  terms  are  collected 


to  obtain 
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0 ■ /N/IH",.,  * Nxy,y  - ”\tt  - R*x,tt)Su° 


+ <Nxy,x  + Ny,y  ' "v!tt  - 


+(Q  + Q - mw  , , + q ) Sw 

x,x  y,y  ,tt  Mz' 


+ (M  + M - Q - I*  a.*.  - Ru° 

x,x  xy,y  x vx,tt  ,tt  yx 


+ (M  + M - Q - I*  . . - Rv  ..)«*  ]dA 
xy » x y,y  *y  *y,tt  ,tt'  V 


♦ H[(Ni  - V"x  + (Niy  * Nxy)ny]5u° 


♦ [ { N ’ - N )n  + (N1  - N )n  ]6v 

xy  xy  x y y y 


+ [ ( Q ’ - Q )n  + ( Q ’ - Q )n  ] <$w 
x x x y y yJ 


+ [ (M • - M )n  + ( M » - M )n  ]5* 
x x x xy  xy  y J yx 


4.20 


+ [ ( M ’ - M )n  + (M1  - M )n  ]6<j»  }ds]dt 

xy  xy  x y y y y' 


Equation  4.20  is  the  variational  form  of  the  plate 
equations  of  motion.  Since  the  6u,  6v,  6w,  5$x,  and 
are  arbitrary  in  the  area,  then  the  coefficients  must 
vanish  independently.  Therefore  the  five  equations  of 


motion  are 
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N + N = mu  , , + Riji 
x,x  xy,y  ,tt  x,tt 


Nxy,x  + Ny,y  = mv,tt  + R*y,tt 


Q +0  + q = mw 

'X|X  y , y Hz  ,tt 


4.21 


M + M - Q ■ I*  4.4.  + Ru  , 
x,x  xy,y  ^x  x,tt  ,tt 


M + M - Q 

xy,x  y,y  ^y 


1^ 


y,  tt 


Equations  4.21  are  the  differential  form  of  the  plate 
equations  of  motion. 

The  boundary  conditions  are  derived  from  the  integral 
around  the  boundary  of  the  element.  On  the  portion  of  the 
boundary  where  the  stress  resultants  are  specified,  the 
variations  are  arbitrary  and  the  coefficients  must  vanish 
to  give  the  natural  boundary  conditions.  On  the  two  sides 
of  a rectangle  where 

n = + 1 n = 0 4.22 

A - y 


the  natural  boundary  conditions  are 


N = N' 
x x 


N = N» 
xy  xy 
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Q = Q * 
x 


4.23 


M = M' 
x x 


M = M' 
xy  xy 


On  the  two  sides  of  a rectangle  where 


nx  . 0 


n = + 1 

y 


4.24 


the  natural  boundary  conditions  are 


N = N' 
xy  xy 


N = N * 

y y 


Q = Q« 

y y 


4.25 


M = M» 
xy  xy 


M = M* 

y y 


On  the  portion  of  the  boundary  where  the  displacements  are 
specified,  5u°,  <5v°,  Sw,  5 x and  are  zero,  giving  the 

essential  boundary  conditions 


u = u ' 


v = v ' 


= w ' 


w 


4.26 
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Comparison  of  the  Variational  and 
Differential  Forms  of  the  Equations 

The  equations  of  motion  of  a plate  have  been  derived 
in  two  forms,  a variational  formulation,  Equation  4.20,  and 
a differential  formulation,  Equation  4.21.  The 
differential  equations  derived  by  variational  methods, 
Equations  4.21,  are  the  same  as  the  differential  Equations 
2.21,  2.22,  2.25,  2.26  and  2.27  derived  from  the  three- 
dimensional  equations  of  motion  in  Chapter  Two.  Since  it 
is  known  that  the  variational  and  differential  formulations 
are  equivalent,  the  two  methods  of  derivation  serve  as  a 
check  on  the  algebra. 

In  the  differential  formulation,  the  problem  is  to 
integrate  a coupled  set  of  partial  differential  equations 
subject  to  known  boundary  conditions.  In  the  variational 
formulation,  the  problem  is  to  find  an  unknown  function  or 
functional  which  minimizes,  maximizes,  or  makes  stationary 
a function  or  system  of  functionals  subject  to  the  same 
given  boundary  conditions.  The  two  problem  formulations 
are  equivalent  because  the  functions  which  satisfy 
differential  equations  and  boundary  conditions  also 
extremize  or  make  stationary  the  functionals  subject  to  the 
same  boundary  conditions. 


CHAPTER  FIVE 

FINITE  ELEMENT  PLATE  EQUATIONS  OF  MOTION 
DERIVED  BY  A VARIATIONAL  METHOD 


Introduction 

The  finite  element  method  is  a procedure  for 
transforming  a problem  from  a variational  or  differential 
form  to  a set  of  linear  algebraic  equations  which  can  be 
solved  by  a computer.  The  finite  element  method  was  not  a 
practical  engineering  method  until  the  high  speed  digital 
computer  was  developed  in  the  1950s.  The  use  of  finite 
elements  has  increased  as  the  cost  of  computer  computations 
has  decreased. 

There  are  many  books  on  the  subject  of  finite 
elements,  but  the  books  become  obsolete  rather  quickly 
because  much  research  is  being  done  in  finite  elements. 

The  subjects  of  plates  and  dynamics  are  usually  covered 
briefly  or  omitted  from  elementary  finite  element  books, 
but  elementary  books  are  necessary  for  someone  beginning 
the  study  of  finite  elemen-ts.  The  following  authors  have 
written  elementary  books  that  provide  a good  general 
background  in  finite  elements:  Desai  (1979),  Hinton  and 

Owen  (1979),  and  Irons  and  Ahmad  (1980),  Becker,  Carey,  and 
Oden  (1981).  The  book  by  Hinton  and  Owen  (1979)  is 
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especially  useful  because  it  has  detailed  information  on 
how  to  write  a finite  element  program.  The  book  by  Huebner 
and  Thornton  (1982)  is  advanced  and  well-written  but  it  is 
oriented  more  towards  fluids  and  heat  transfer. 

The  following  books  cover  plates  and  dynamics  in 
detail.  The  book  by  Zienkiewicz  (1977)  is  the  most 
comprehensive  on  plates  and  dynamics  but  is  challenging  to 
read.  Other  books  which  cover  plates  and  dynamics  are  Cook 
(1981)  and  Reddy  (1984).  The  most  complete  coverage  of 
dynamics  and  solution  of  dynamic  equations  is  the  book  by 
Bathe  (1982). 


Dynamic  Finite  Element  Equations 
The  finite  element  equations  of  motion  for  a laminated 
plate  will  be  derived  in  this  chapter  by  using  Hamilton's 
principle.  The  variational  statement  of  the  strain  energy, 
the  work,  and  the  kinetic  energy  from  Chapter  Four  will  be 
used. 

The  finite  element  equations  of  motion  for  an  element 
have  the  form 


[ k ] { d } + [m]{d^tt|  = { f } 5.1 

where  [k]  is  the  element  stiffness  matrix,  [m]  is  the 
element  mass  matrix,  fd}  is  the  element  nodal  displacement 


and  rotation  vector,  {d  is  the  element  nodal 
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acceleration  vector  and  {f}  is  the  element  nodal  force 
vector.  The  set  of  linear  equations,  Equation  5.1,  is 
solved  for  the  unknown  displacements  { d } and  accelerations 
{dftt}‘  Since  all  the  unknowns  are  displacements  or 
rotations,  this  is  the  displacement  finite  element 
formulation . 

The  displacement  vector  contains  the  five  element 
degrees  of  freedom, 

{u} 

M 

{ W } 

} 

1 X 

{*  } 
y 

Each  vector  contains  the  displacement  at  each  of  the  eight 
nodes  (for  an  eight  mode  element). 

T 

{u}  = [u1  u2  u5  u4  u5  u6  Uy  Ug]  5.3 

The  other  vectors  { v } , {w},  and  { <i»y } are  similar. 

The  nodal  force  vector  contains  the  applied  forces  and 


moments , 


55 


Each  subvector  {f}  contains  one  force  or  moment  acting  at 
each  node  point. 

The  stiffness  matrix  contains  submatrices, 


[k11]  [k12] 

[ 0 ] 

tk1*] 

[k15] 

[k22] 

[ o ] 

tk2*] 

tk25] 

[k55] 

tk34) 

!k35] 

5.5 

symmetric 

tk44] 

[k45] 

[k55] 

Each  submatrix  is  square  with  the  number  of  rows  and 
columns  equal  to  the  number  of  node  points  in  the  element 
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The  element  mass  matrix  has  the  form 


Cm  ] 0 


0 


[m]  = 


0 

0 

0 

0 


[m  ] 0 


0 


[m1]  0 


0 


[in  ] 0 


0 [rn  ] 


5.6 


where  the  submatrix  [m"'  ] contains  the  mass  terms  and  the 
submatrix  [m^]  contains  the  rotational  moment  of  inertia 
terms.  The  mass  matrix  and  subraatrices  are  the  same  size 
as  the  stiffness  matrix  and  submatrices. 

The  finite  element  stiffness  and  mass  matrix  will  be 
derived  in  this  chapter  for  the  symmetric  layup  cross-ply 
laminate  discussed  in  Chapter  Two.  For  a cross-ply 
laminate 


'16 


= D 


26 


= 0 


For  a symmetric  layup  laminate 


and 


5.7 


5.8 


R = 0 


5.9 
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The  Bj^  terms  couple  the  inplane  and  out-of-plane 
displacements.  Since  B^j  is  zero  for  a symmetric  laminate, 
the  two  inplane  equations  are  uncoupled  from  the  three  out- 
of-plane  equations.  If  there  are  no  inplane  loads,  then 
the  are  not  needed  to  solve  the  problem.  For  a cross- 
ply  symmetric  layup  laminate,  the  force  displacement 
relationships,  Equations  2.18,  reduce  to 


M = D.,.^  + D.0<|> 

x 11  x,x  12  y,y 


y 1 2 x , x 22  y,y 


M = Drc  ( »l>  + ^ ) 

xy  66  y,x  x,y 


Q = AC(.(w  + ip  ) 
x 55  ,x  rx 


Q = A . . (w  + t|/  ) 

y 44  ,y  y 


5.10a 
5.10b 
5.10c 
5. 10d 
5. 10e 


Finite  Element  Stiffness  Matrix 
The  element  stiffness  matrix  will  be  derived  from  the 
variational  statement  of  the  strain  energy,  Equation  4.7, 
by  a method  outlined  by  Reddy  (1980).  The  force- 
displacement  relations,  Equations  5.10,  are  used  to  obtain 


58 


6U 

= //[(D^x.x  + 

°1 2^ 

y.y 

)5*x,x 

+ 

(D12*x,x  * 

D22^y 

,y) 

6 ip 

y.y 

+ 

D66^y,x  * 

*x,y^ 

6 

( d;  + di  ) 

y > x *x,y' 

5 

.11 

+ 

A55(w,x  + * 

*>  5 

(w 

t 

+ <p  ) 
x yx 

+ 

A44(\y  * * 

y}  5 

(w 

* 

y + 4<y ) ]dxdy 

basis 

of 

the  finite 

element 

method  is  to 

approximate 

unknowns 

w,  i|»x,  and 

*y  as 

a 

sura 

n 

w 

(x,y,t)  = z 
d 

\]<x' 

y) 

w .(t) 
J 

5. 

12a 

n 

x(x,y,t)  = z 
3 

4j(X 

,y) 

*xj(t) 

5. 

12b 

n 

(x,y,t)  = E 

*i(x 

»y ) 

*ya(t) 

5. 

12c 

where  $j(x,y)  are  shape  functions  and  Wj(t),  j ( t ) , and 
j ( t ) are  unknown  displacements  and  rotations  at  each  node 
of  the  elementjand  n is  the  number  of  nodes  in  the 
element.  The  shape  functions  have  a value  of  one  at  the 
jth  node  and  zero  at  the  other  nodes.  The  summation  is 
over  the  number  of  nodes  in  the  element.  The  partial 
derivatives  are 
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n 

w = £ <t>  . w , 
>x  ^ 0.x  J 


n 

w = £ $ . w . 

»y  . o»y  j 


n 


4>  = £ $ . di 

x»x  J j.X^XJ 


*x>y  ?j  ^j.y^xj 


n 

ip  = E $ , il; 
y»x  j J,xvyj 


n 

l|)  = £ $ d/ 

y.y  j J.y^yJ 


The  variations  are 


n 

6w  = £ 


6*X  = l V*Xi 


n 

64  = £ . 

y i i yi 


The  following  integrals  are  defined 


S°3  = /#  i*jdxdy 


Sij  * '‘l,x‘;),ydltdy 


Sij  - ^t,x‘jdXdy 


SiJ  = /*l.y*ddxdy 


Sij  * /Jl,x*j,xdxdy 


SH  = /♦ . t . dxdy 

ij  J i»y  o,y  y 


The  variation  of  the  strain  energy,  Equation  4.7, 


5.13a 


5.13b 


5.13c 


5.14a 


5.14b 


5.14c 


5.15 


with 


the  finite  element  approximations,  Equations  5.12,  5.13, 
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5.14,  and  the  definitions  of  Equations  5.15,  yields 


n n 


6U  = E E {(D  SXX'I»  + D SX^4>  ) 6ii) 
u * . uu11aijvx;J  u12  l^yj'  xi 


+ (D„_S?X<I>  . + d00sH*  -)6*  < 
12  13  xj  22  lj  yj  yi 


5.16 


* D66t'Si?ya  * S«V5V  * (Sid*yd  * Sij*xd)S*xi] 

* * ‘S?h  * 

* *44“sl$"i  * sij*y3)6wi  * (sid“j  * S“W» 


In  matrix  form 


«U  = [[«w],[«*x],[6*y]] 


[k43]  [k44]  [k 


5.17 


The  submatrices  of  the  stiffness  matrix  are 


[k33]  . 455[Sxx]  + A44[Syy] 
[k34]  - A55[SXO] 

[k35]  - A44[Sy°] 


5.18 
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[k44] 

= D1 1 [SXX  ] 

+ d66[s") 

+ a55[s°°] 

[k45] 

= D12tsxn 

* D66tSyX] 

[k55] 

r— i 

>> 

CO 

CM 

CM 

O 

II 

+ d66[s-] 

♦ a44[s°°] 

The  stiffness  matrix  can  be  shown  to  be  symmetric. 
Equations  5.15  can  be  summarized  as 


39.  39  . 

ff  jttJ1  dxdy 


5.19 


Since  the  order  of  multiplication  can  be  reversed 


In  matrix  form 


[S5n]  = [Sn5]T 


5.20 


5.21 


Therefore  the  stiffness  matrix  is  symmetric  and 

[k43i  = tk34]1 

[k53]  . [k33]T 


5.22 


62 


Finite  Element  Force  Vector 
The  element  nodal  force  vector  will  be  derived  from 
the  virtual  work  done  by  distributed  loads  on  the  plate, 
Equation  4.12,  and  the  virtual  work  done  by  point  loads  on 
an  element  in  the  plate, 


n 

6w  = - E F. 6w.  5.23 

i 1 1 

The  finite  element  approximations,  Equations  5.12,  are 
substituted  into  the  virtual  work,  Equations  4.12  and  5.23, 
to  yield 


n 


>W  = E - [(Q'n  + Q'n  )$,<5w.  + (M'n  + M'  n . 

J xx  y y i i x x xy  y l xi 


+ ( M * n + M'n  .Ids 

xy  x y y l yi J 


//qz$i6widxdy 


n 

E F . 6w. 
i 1 1 


In  matrix  form 


(f3} 


5.24 


5W  = - [[6w],[«*v],[«<,y]]  < (f4) 

(f5) 


5.25 


where 


63 


= *(Qxnx  + Qyny)*ids  + //qz$idxd  y + Fi  5.26a 

fi  = *(Mxnx  + Mxyny)$ids  5.26b 

f-  = 4>(M;ynx  + M^ny)$.ds  5.26c 


Finite  Element  Mass  Matrix 
The  finite  element  mass  matrix  will  be  derived  from 
the  variational  statement  of  the  kinetic  energy,  Equation 
4.19*  The  partial  derivatives  with  respect  to  time  of  the 
finite  element  approximations,  Equations  5.16,  are 


n 

w , . = E $ . w . , , 
»tt  j J J,tt 


n 

^x , tt  * *xj,tt 

J 


n 

♦y.tt  - £ *j 


5.27 


Equations  5-27  are  substituted  into  the  variational 
statement  of  the  kinetic  energy,  Equation  4.19*  to  obtain 


t.  n n t1 

6 / 1 T dt  = E E / 

t0  1 d fc0 


S°°[IU  • . 

i JL  xj,tt  vxi 


+ 


'i'  - . . 5iJ»  . ) 
yj*tt  yyi' 


+ m w . ,,  6w. ]dt 
O.tt  iJ 


5.28 
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The  rotational  coupling  inertia,  R,  is  zero  for  a symmetric 
layup  plate.  The  negative  sign  is  a result  of  integration 
by  parts  and  cancels  with  the  negative  signs  in  the  strain 
and  the  work  in  Hamilton’s  principle. 

In  matrix  form 


T dt 


[[«w] , [6*XJ , [6*y]][m] 


5.29 


^y , tt  ^ 


where 


[m] 


m[S°°  ] 0 0 

0 I[S°°]  0 

0 0 I [S00 ] 


5.50 


The  mass  per  unit  area,  m,  and  the  rotary  inertia  per  unit 
area,  I,  are  defined  by  Equations  2.23  and  2.38 
respectively. 


Hamilton’s  Principle 

Hamilton's  principle  is  used  to  find  the  dynamic 
finite  element  equations.  The  variational  statements  of 
the  stiffness,  Equation  5.17,  kinetic  energy,  Equation 
5.29,  and  work,  Equation  5.25,  are  substituted  into 
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Hamilton's  principle,  Equation  4.1,  to  obtain 


/ 


t 


0 


[[«*], C«*x],[«*y]] 


= 0 


5.31 


Since  this  relation  is  true  for  every  virtual  displacement 
then 


Equation  5.32  is  the  finite  element  equation  for  motion  for 
a single  element.  The  stiffness  matrix  [k]  is  given  in 
Equation  5.18,  the  force  vector  [f]  is  given  in  Equation 
5. 26, and  the  mass  matrix  [m]  is  given  in  Equation  5.30. 

The  displacements  {w} , { ^ x } , and  {♦y}  are  the  unknowns.  A 
numerical  method  presented  in  Chapter  Seven  is  used  to 
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integrate  the  accelerations  {w  tt} , {’J»x  tt},  and  {♦y  ttl 
with  respect  to  time. 


CHAPTER  SIX 

FINITE  ELEMENT  PLATE  EQUATIONS  OF  MOTION 
DERIVED  BY  GALERKIN  * S METHOD 

Introduction 

The  equations  of  motion  of  a plate  can  be  expressed  in 
either  a differential  formulation  or  a variational 
formulation,  which  are  equivalent.  The  variational  form  of 
the  plate  equations  of  motion  was  derived  from  Hamilton's 
principle  in  Chapter  Four  and  used  in  Chapter  Five  to 
derive  the  finite  element  plate  equations  of  motion.  The 
differential  form  of  the  plate  equations  of  motion  was 
derived  from  the  three-dimensional  equations  of  motion  of 
an  elastic  body  in  Chapter  Two  and  from  the  variational 
form  in  Chapter  Four.  In  this  chapter  a few  terms  of  the 
finite  element  plate  equations  of  motion  will  be  derived 
from  the  differential  form  of  the  plate  equations  of  motion 
by  Galerkin's  method.  The  purpose  of  this  chapter  is  to 
derive  the  finite  element  plate  equations  of  motion  by  a 
second  method  which  provides  a check  on  the  results  derived 
in  Chapter  Five. 
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Plate  Equations  of  Motion  in  Terms  of  Displacement 
The  differential  form  of  the  plate  equations  of  motion 
of  a plate  are  stated  in  Equation  4.21  in  terms  of  stress 
resultants,  Mv,  Mv__,  Qv,  and  Qw.  The  force- 

a y Ay  a y 

displacement  relations  for  a cross-ply  symmetric  layup 
laminate,  Equations  5.10,  are  substituted  into  Equation 
4.21  to  obtain  the  plate  equations  of  motion  in  terms  of 
displacements.  The  results  can  be  put  in  the  following 
form 


[L]  {d}  + {f}  = 0 


6.1 


where  L is  a matrix  of  differential  operators 


The  individual  differential  operators  are 


6.2 


L33 

- a55( 

' , XX 

+ A44(  >,yy  - m( 

^ ,tt 

L34 

* A55< 

>.« 

L43  * 

-A55< 

L35 

' A44< 

\x 

L53  - 

-A44( 

L44 

' -S55( 

) * 

D11 ( >,xx*D66< 

>yy 

6.3 
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Equation  6.3  is  a special  case  of  the  more  general  plate 
equations  of  motion  in  term  of  displacement  stated  by  Yang, 
Norris,  and  Stavsky  (1966,  p.  670).  Note  that  the 
differential  operators,  Equations  6.3,  have  the  same 


components  of  the  finite  element  stiffness  matrix,  Equation 
5.18. 

The  finite  element  plate  equations  will  be  derived 
from  Equations  6.2  and  6.3  by  the  Galerkin  method,  which  is 
a particular  type  of  the  method  of  weighted  residuals.  The 
Galerkin  method  is  equivalent  to  the  variational  method 
used  in  Chapter  Five.  The  weak  form  of  the  equilibrium 
equations  discussed  by  Zienkiewicz  (1977,  p.  45)  is  also 
equivalent  to  the  variational  and  Galerkin  methods. 

The  finite  element  equations  are  derived  by 
substituting  the  finite  element  approximations,  Equations 
5.12,  into  the  operator  Equations  6.1  and  6.2,  which 
results  in  a residual  (an  error) 


The  error  is  a function  of  the  independent  variables  and 


material  properties,  and  as  the  corresponding 


6.4 
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to  n,  the  number  of  nodes  in  the  element.  In  the  weighted 
residual  method,  the  parameters  are  found  by  multiplying 
the  error  by  a weight  function,  integrating  over  the  area, 
and  setting  the  integral  equal  zero. 


//♦i(x,y)(L(x,y,t,w^,^xj,*yj)  + f^)dxdy  = 0 


6.5 


Since  the  same  functions  $ are  used  for  both  the  assumed 
solution  and  the  weight  function,  this  method  is  known  as 
the  Galerkin  method. 


.35 


Three  terms  in  the  stiffness  matrix,  k^,  ^ij*  an^ 


k. .,  will  be  derived  from  the  differential  operator  terms, 
^33’  ^34*  an^  ^35’  illustrate  the  Galerkin  method.  The 
terms,  L33,  ^34*  and  L35,  defined  in  Equations  6.2  and  6.3 
are  substituted  into  Equation  6.5  to  obtain 


//$.[Accw  + A..w  + Acci|»  + A. 

iL  55  , x x 44  ,yy  55  x,x  44  y,y 


- m w , . + q ]dxdy  = 0 

y U o Z 


6 . 6 


Integration  by  parts  is  performed  using  Equation  4.9  and 


ffF  dxdy  = - JJ  -§|  -§|  dxdy  + fF  nxds  6.7a 

3x 

2 

t 3 G j j f r 3F  3G  , , r,-,  3G 

// F ^2  dxdy  = " " 17  T7  dxdy  + fF  17  nyds 


6.7b 
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where  n 

normal 

is 


and  n 
x y 

and  the 


are  the  direction 
x and  y directions 


cosines  between  the 
respectively.  The  result 


//[-A55*t,x\x  - A44*i,y\y  ' A55*l,x*x  ' A44»l,y*y 

- * *ldzldxd!''  * #»ltA55*,*nx  * A44l',yny  6'3 

+ Act-ih  n + A. .ib  n ]ds  = 0 
55vx  x 44vy  y 

The  finite  element  approximations,  Equations  5.12  and  5.13, 
and  the  definitions  of  integrals.  Equations  5.15,  are  used 
to  obtain 


A55SiJwj  + A44Sijwj  + A55^i j^x j + A44Sij*yd  + mSijwd,tt 
- //*i^zdxdy  + K[A55(w,x  + *x)nx  + A44 (w,y  + Vny]ds 


6.9 


If  Equation  6.9  is  put  in  the  matrix  equation  form  defined 
by  Equations  5.1,  5.2,  5.3,  5.4,  and  5.5,  the  stiffness 
matrices  are 


[k33]  = a55[sxx]  * A44[Syy] 

[k34]  - A55[Sx°] 

[k33]  . a44[s»»] 


6.10 
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which  is  part  of  Equation  5.18.  The  mass  matrix  is 


6.11 


which  is  part  of  Equation  5.30.  Equations  5.10  d and  e are 
applied  to  Equation  6.9,  and  the  results  are  put  into  the 
form  defined  by  Equations  5.1,  and  5.4.  The  point  loads, 
F^,  are  added  to  obtain  the  components  of  the  force  vector 


which  is  Equation  5.26  a. 

The  other  terms  in  the  stiffness  and  mass  matrices  and 
force  vector  are  obtained  by  applying  Galerkin's  method 
with  the  other  terms  in  the  differential  operator.  The 
same  method  is  used  and  the  results  agree  with  the  results 
in  Chapter  Five,  but  the  algebra  will  not  be  shown  here. 

In  this  chapter,  it  has  been  shown  that  the 
variational  and  differential  form  of  the  plate  equations  of 
motion  yield  the  same  finite  element  equations.  Less 
algebra  is  required  with  the  variational  form  of  the 
equations  because  integration  by  parts  is  not  required. 

For  some  types  of  problems  a variational  form  is  unknown  or 
nonexistent,  and  the  differential  form  must  be  used. 


f . 

l 


= f$.(Q’n  + Q’n  )ds  + //$. q dxdy  + F. 

■>  i x x y y i z J x 


6.12 


CHAPTER  SEVEN 

DYNAMIC  SOLUTION  OF  THE  FINITE  ELEMENT  EQUATIONS 


Newmark's  Method  for  Integration  with  Respect  to  Time 
The  finite  element  equations  of  motion  for  a laminated 
plate  were  derived  from  Hamilton's  principle  in  Chapter 
Five.  The  element  matrices  and  vectors  in  the  equations  of 
motion  for  one  element,  Equation  5.1,  are  assembled  into 
the  global  matrices  and  vectors  for  the  entire  plate.  The 
global  equations  of  motion  have  the  form 

[K]{D}  + [M]{D>tt}  = {F}  7.1 

where  [K]  is  the  global  stiffness  matrix,  [M]  is  the  global 
mass  and  moment  of  inertia  matrix,  [D]  is  the  global  nodal 
displacemnt  and  rotation  vector,  [D^^.]  is  the  global  nodal 
acceleration  vector,  and  {F}  is  the  global  nodal  force 
vector.  The  stiffness  matrix  [ K ] contains  functions  of  the 
shape  function,  $,  which  depend  only  on  spatial 
coordinates.  Equation  7.1  is  a set  of  ordinary  linear 
differential  equations  and  was  found  by  making  a spatial 
approximation.  Such  a procedure  is  called  semidiscrete 
approximation . 
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Time  discretization  was  performed  by  Newmark's 
method.  The  displacement  and  acceleration  vectors  {D}  and 
{D  are  functions  of  time  and  are  found  by  integrating 
Equation  7.1  numerically  with  respect  to  time. 

Newmark's  equation,  according  to  Cook  (1981,  p.  322), 
has  the  form 

[[K]  ♦ a0[M]J(D}n+1  - (F}n+1  * [M]{a0(D}n  * a1)Dj(.}n 


+ a2<D,ttU 

7.2 

= 1/B(At)2 

7.3a 

= i/e(At) 

7.3b 

= 1/28  - 1 

7.3c 

The  parameters  a and  S control  the  stability.  For 

a = 1/2  B = 1/4  7.4 

Newmark's  method  is  an  implicit  method  and  is  uncondi- 
tionally stable  for  any  time  step  and  introduces  no 
artificial  damping  according  to  Cook  (1981,  p.  322).  The 
displacement  at  time  step  n+1  is 


75 


{°>n*1  * f°n>  * )2t(i-  S){Dttt)n 

+ S(D,tt}n+11  7'5 


The  velocity  at  time  step  n+1  is 


<D,t>n+1  - (D,tln  * t «(0>tt}n+1]  7.6 


The  acceleration  at  time  step  n+1  is 


(D  , lDUl-<Dln-<atHD.t)n 
1 6(At)2  + 6 


7.7 


The  algorithm  operates  as  follows.  At  n = 0 the 
initial  conditions  prescribe  { D } q , {D  ^.}q,  and  {D  ^Jq. 
Equations  7.2  and  7.3  are  solved  for  the  displacements 
{D } 1 . Equations  7.6  and  7.7  are  solved  for  velocity  and 
accelerations  {D^}^  and  The  process  is  repeated 

for  n = 1.  Equations  7.2  and  7.3  are  solved  for  the 
displacement  { D } 2 • Then  Equations  7.6  and  7.7  are  solved 
for  {D  j.}2  and  { D ^ } 2 • The  process  is  continued  until  the 
displacements  at  the  desired  time  are  found. 

Newmark's  method  with  values  of  alpha  and  beta  used 
here  is  stable  for  all  time  steps.  Therefore,  there  is  not 
a critical  minimum  time  step  required  for  stability. 
Different  time  steps  were  tried  and  results  showed  that  the 
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displacements  and  stresses  converge  towards  a particular 
value  as  the  time  step  is  reduced. 

Mass  Matrix 

The  element  mass  matrix  in  Equation  5-30  is  a full 
matrix  and  is  called  a consistent  mass  matrix.  There  are 
two  other  types  of  mass  matrices  called  the  diagonal  and 
lumped-mass  matrix  which  have  non-zero  terms  only  on  the 
main  diagonal.  The  lumped  mass  matrix  is  formed  by  lumping 
the  mass  equally  at  each  node.  The  diagonal  mass  matrix  is 
formed  by  taking  the  diagonal  terms  of  the  consistent  mass 
matrix  and  scaling  the  terms  so  that  the  total  mass  of  the 
element  is  preserved. 

The  choice  of  consistent,  diagonal,  or  lumped  mass 
depends  on  the  particular  problem  being  solved.  For  a 
simply  supported  thick  square  plate,  the  diagonal  mass 
matrix  gives  more  accurate  natural  frequencies  than  either 
the  consistent  or  the  lumped  mass  matrix  for  the  majority 
of  the  modes  according  to  Hinton,  Rock,  and  Zienkiewicz 
(1976).  Therefore,  the  diagonal  mass  matrix  was  used. 

The  finite  element  mass  matrix,  Equation  5*30,  is  a 
consistent  mass  matrix,  which  was  replaced  by  a diagonal 
mass  matrix.  The  total  mass  of  the  element,  rag , 
to  Equation  2.24,  is 


according 
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m 


e 


p h A 


7.8 


where  p is  the  density,  mass  per  unit  volume,  h is  the 
thickness  of  the  plate,  and  A is  the  area  of  the  element. 
For  the  diagonal  mass  matrix,  this  total  mass  is 
distributed  among  the  nodes  of  the  eight-node  element  as 
follows.  The  mass  assigned  to  each  corner  node  is 


as  shown  in  Figure  7.1. 

The  total  moment  of  inertia  of  the  element,  according 
to  Equation  2.29,  is 


7.9a 


and  the  mass  assigned  to  each  side  node  is 


7.9b 


I 


e 


7.10 


This  is  distributed  as  follows.  The  moment  of  inertia 
assigned  to  each  corner  node  is 
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p h5  A 
c " 432 


7.11a 


The  moment  of  inertia  assigned  to  each  side  node  is 


p h^  A 
s 54 


7.11b 


When  the  diagonal  masses  are  assembled  into  the  global  mass 
matrix,  the  total  mass  at  each  side  node  of  the  element  is 
4 times  the  total  mass  at  each  corner  node. 


Impact 

The  purpose  of  this  study  was  to  understand  the  impact 
response  of  a plate.  Two  types  of  impact  were  considered, 
an  inelastic  impact  and  an  elastic  impact.  For  the 
inelastic  impact  the  projectile  hits  the  plate  and 
sticks.  It  does  not  penetrate  the  plate  or  bounce  off. 

For  the  elastic  impact  the  projectile  hits  the  plate 
surface,  deforms  the  plate  surface,  and  bounces  off.  The 
force  between  the  projectile  and  the  plate  is  governed  by 
the  deformation  interaction  between  the  projectile  and 
plate.  The  actual  interaction  law  is  unknown,  but  for  this 
study  a Hertzian  elastic  contact  law  was  assumed.  The 
results  of  the  two  types  of  impact  will  be  compared  to 
determine  the  effect  of  the  type  of  impact  on  the  plate 
displacement  and  stresses. 
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Figure  7.1 


Distribution  of  the  Element  Mass  and  Moment  of 
Inertia  among  the  Eight  Nodes  for  a Diagonal 
Mass  Matrix 
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Inelastic  Impact 

In  this  study  an  inelastic  impact  is  modeled  as  an 
inelastic  collision  between  the  projectile  and  the  center 
finite  element  node  with  a lumped  mass  as  suggested  by 
Goldsmith  (I960,  p.  56).  The  plate  is  initially  at  rest. 
After  the  impact  the  center  node  of  the  plate  and  the 
projectile  move  at  the  same  velocity.  Momentum  is 
conserved,  which  means  that  the  momentum  of  the  projectile 
before  impact  is  equal  to  the  momentum  of  the  projectile 
and  the  mass  of  the  impacted  node  after  impact. 


m v 
P P 


m ) 

n 


0 


7.12 


where  mp  is  the  mass  of  the  projectile,  mn  is  the  mass  of 
the  node  where  impact  occurs,  Vp  is  the  velocity  of  the 
projectile  before  impact,  and  vq  is  the  velocity  of  the 
projectile  and  the  impacted  node  after  impact.  Equation 
7.12  can  be  solved  for  Vq. 


v 


P 


V 


7.13 


For  a given  mass  and  velocity  of  the  projectile,  mp  and  Vp, 

the  computer  program  calculates  Vq,  which  is  used  as  the 

initial  velocity  of  the  center  node  of  the  plate  in  the 

finite  element  program.  The  computer  adds  the  mass  of  the 

projectile,  m , to  the  mass  of  the  center  node  and  uses 
P 
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this  as  the  mass  of  the  center  node  in  the  subsequent 
calculation  of  the  motion  of  the  plate. 

The  mass  of  the  projectile,  mp,  is  much  greater  than 
the  mass  of  the  node,  mn,  which  models  the  mass  of  the  area 
of  impact  of  the  plate.  Since  mp  is  much  greater  than  mn, 
the  change  in  velocity  of  the  projectile  before  and  after 
impact  is  small  according  to  Equation  7.13,  which  means  the 
energy  loss  is  small.  Most  of  the  kinetic  energy  of  the 
projectile  is  transferred  to  the  plate  and  only  a small 
amount  of  energy  is  lost  at  impact. 

Elastic  Impact 

The  elastic  impact  included  indentation  of  the  plate 
surface  by  the  projectile.  The  force  between  the 
projectile  and  the  plate  during  impact  was  assumed  to  be 
governed  by  a law  for  contact  between  two  bodies  of  the 
form 


F = H (r  - w )p  7.14 

where  F is  the  force,  H is  a material  constant,  r is  the 
displacement  of  the  projectile  and  w is  the  displacement  of 
the  center  of  the  plate,  and  p is  an  exponent.  Both 
displacements  are  measured  from  the  surface  of  the 
undeformed  plate. 
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The  deformation  interaction  between  the  projectile  and 
the  plate  in  the  actual  interaction  is  unknown,  but  for 
this  study  a Hertzian  elastic  contact  with  p = 1.5  was 
assumed.  Hertzian  theory  of  impact  is  discussed  by 
Goldsmith  (I960,  p.  37)  and  Love  (1944,  p.  199).  The 
static  force  required  to  indent  glass-epoxy  and  graphite- 
epoxy  laminates  with  a round  steel  ball  was  determined 
experimentally  by  Sun,  Sankar,  and  Tan  (1981).  They  found 
that  a Hertzian  elastic  contact  law  (p  = 1.5)  best  fit  the 
loading  curve.  The  Hertzian  contact  law  was  used  in  this 
study  for  dynamic  loading  and  unloading  between  a cylinder 
and  a glass-epoxy  plate. 

The  position,  r,  and  the  velocity,  r +.  of  the 

9 ^ 

projectile  are  calculated  using  laws  of  dynamics  of  a 
particle  under  a constant  force  expressed  as 


r , = r ^ + — At 

n + 1 , t n,  t nip 


7.15 


rn+1  * rn  + rn,t  St  + 25“  <4t>2 

P 


7.16 


where  mp  is  the  mass  of  the  projectile,  At  is  the  time 
step,  and  F is  the  force  acting  on  the  projectile  caused  by 
elastic  impact.  The  subscripts  n and  n+1  indicate  time 
steps.  At  n = 0,  the  initial  position,  rQ,  is  zero  because 
the  projectile  is  initially  touching  the  surface  of  the 
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undeformed  plate.  The  initial  velocity,  rQ  ^ , is  known  and 
the  initial  force,  FQ,  is  zero. 

At  the  end  of  the  first  time  step,  the  projectile 
position,  r ^ , and  velocity,  r^^.,  are  calculated  using 
Equations  7.15  and  7.16.  The  projectile  has  indented  the 
surface  of  the  plate,  and  Equation  7.14  is  used  to 
calculate  the  force,  F1 , between  the  projectile  and  the 
plate.  The  force  is  used  in  the  finite  element  program  to 
calculate  the  plate  deflection  during  the  next  time  step. 

At  n equals  1,  Equations  7.15  and  7.16  are  used  to 
calculate  the  projectile  position,  r ?,  and  velocity, 
r2,t*  This  process  is  continued  for  each  time  step. 


CHAPTER  EIGHT 

STEPS  IN  A FINITE  ELEMENT  ANALYSIS 


A typical  dynamic  finite  element  analysis  of  an 
elastic  body  consists  of  the  following  steps:  (1) 

Discretization  of  the  domain  into  finite  elements,  (2) 
Selection  of  element  shape  functions,  (3)  Formation  of  the 
element  stiffness  and  mass  matrices,  (4)  Assembly  of  the 
element  matrices,  (5)  Imposition  of  the  boundary  conditions 
and  static  or  dynamic  forces,  (6)  Solution  of  the  dynamic 
finite  element  equations,  and  (7)  Calculation  of  the 
stresses  and  strains. 

The  first  step  is  discretization  of  the  plate  into 
subregions  called  finite  elements.  The  finite  element 
program  developed  in  this  study  solves  the  problem  of  a 
square  plate  impacted  at  the  center.  Square  elements  were 
chosen.  The  locations  of  the  elements  for  examples  using 
1 » 4,  9,  25,  36,  and  64  elements  are  shown  in  Figure  3.1. 
Since  the  plate  geometry,  boundary  conditions,  materials, 
and  loading  are  symmetric  with  respect  to  axes  through  the 
center  of  the  plate,  and  the  layup  is  symmetric  with 
respect  to  the  midplane,  it  is  sufficient  to  analyze  only 
one  quarter  of  the  plate. 
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(a) 

(c) 

(e) 


(d) 


Figure  8.1.  Placement  of  Elements  on  the  Plate 
a)  One  Element,  b)  Four  Elements, 
c)  Nine  Elements,  d)  Twenty-five  Elements, 
e)  Thirty-six  Elements,  f)  Sixty-four 
Elements 
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Step  two  is  selection  of  element  shape  functions.  A 
quadratic,  serendipity  element  was  chosen.  The  element  is 
a square  (or  rectangular),  two  dimensional  eight  node 
element  as  shown  in  Figure  8.2.  There  are  three  unknowns 
at  each  node  point:  displacement,  w,  and  rotations,  <t>x  and 

fy.  Each  element  has  24  unknowns.  Therefore,  the  size  of 
the  element  stiffness  matrix  is  24  x 24. 

There  are  eight  quadratic  shape  functions,  one 
assigned  to  each  node. 


^ = 1/4  ( 1 -x) ( 1 -y ) ( -x-y-1 ) 
*2  = 1/4  ( 1 +x ) ( 1 -y ) ( +x-y-1 ) 
$3  = 1/4  (1 +x) (1+y) (+x+y-1 ) 
% = ( 1 -x ) ( 1 +y ) ( -x+y-1 ) 

*5  = 1/2  (l-x2)(1-y) 

*6  = 1/2  (1 +x ) (1 -y2 ) 

*7  = 1/2  ( 1 -x2 ) ( 1 +y ) 

*8  = 1/2  (1-x)(1-y2) 


8.1 


Each  shape  function  has  the  value  one  at  its  node  and  is 
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zero  at  all  the  other  nodes.  One  of  the  shape  functions  is 
shown  in  Figure  8.3.  The  value  of  the  shape  function  is 
represented  by  the  height  perpendicular  to  the  x-y  plane. 

A quadratic  shape  function  fits  a curved  surface 
through  the  node  points  as  shown  in  Figure  8.3.  Linear 
shape  functions  that  fit  straight  lines  through  the  node 
points  are  simpler  to  program  but  were  not  used  for  an 
important  reason:  for  the  same  number  of  node  points  in 

the  global  equations,  a linear  shape  function  may  produce 
displacements  that  are  as  accurate  as  displacements 
produced  by  quadratic  shape  functions,  but  the  stresses  and 
strains  will  not  be  as  accurate  because  these  depend  on  the 
derivatives  of  the  shape  function  with  respect  to  x and  y. 

The  third  step  is  formation  of  the  element  stiffness 
and  mass  matrices  and  the  force  vector.  The  stiffness 
matrix  is  formed  from  Equation  5.18.  The  element  mass 
matrix  is  formed  from  Equations  5.30,  7.9,  and  7.11,  and 
the  element  force  vector  is  formed  from  Equations  5.26. 

The  integrals  in  Equation  5.15  were  found  by  numerical 
integration  using  Gauss  Quadrature. 

The  fourth  step  is  assembly  of  the  element  stiffness 
matrix,  [k],  the  diagonal  mass  matrix,  [m],  and  the  force 
vector,  { f } , into  the  global  stiffness  matrix,  [K],  global 
mass  matrix,  [M],  and  global  force  vector,  {f}. 

The  fifth  step  is  imposition  of  boundary  conditions. 
The  boundary  conditions  for  a simply  supported  plate  are 
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Figure  8.2.  A Quadratic  (8  node)  Serendipity  Element 


Figure  8.3.  A Typical  Quadratic  Shape  Function 
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shown  in  Figure  8.4a  and  for  a clamped  plate  are  shown  in 
Figure  8.4b.  The  boundary  conditions  along  the  axes  are 
chosen  to  produce  displacements  which  are  symmetric  with 
respect  to  the  axes. 

The  sixth  step  is  solution  of  the  set  of  linear 
equations,  Equation  7.1,  which  are  the  result  of  the  finite 
element  process.  The  equations  are  solved  by  Gauss 
elimination  and  Newmark's  method  presented  in  Chapter 
Seven.  The  results  are  displacement,  w,  and  rotations,  <l<x 
and  ^y,  at  each  node. 

The  seventh  and  last  step  is  calculation  of  stress  and 
strain.  The  strains  are  calculated  from  the  displacements 
by  using  Equation  2.7  and  the  stresses  are  calculated  from 
the  displacements  by  using  Equation  2.13.  The  program 
produces  graphs  of  displacement  and  stress  versus  distance 
from  impact  for  various  times  which  are  shown  in  Chapter 


Ten . 
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Figure  8.4.  Boundary  Conditions  for  a Plate 

a)  Simply  Supported,  b)  Clamped 


CHAPTER  NINE 

VERIFICATION  BY  SOLVING  PROBLEMS 
WITH  KNOWN  SOLUTIONS 

Introduction 

A finite  element  program  is  not  useful  until  it  has 
been  verified  to  give  correct  answers.  The  finite  element  * 
program  developed  in  this  study,  which  will  be  referred  to 
as  Petersen's  program,  was  verified  by  solving  problems 
which  have  analytical  solutions  and  by  comparing  with  well 
established  finite  element  programs.  The  following 
problems  with  analytical  solutions  were  considered:  static 

deflection  of  an  isotropic  plate  with  a point  load  at  the 
center,  cylindrical  bending  of  an  isotropic  plate,  and  free 
vibration  of  an  isotropic  plate.  It  was  shown  that  the 
finite  element  solution  converges  toward  the  correct  answer 
as  the  number  of  elements  increases  from  1 to  25. 

Plate  Material  Properties,  Weights,  and  Dimensions 
The  material  properties,  weights,  and  dimensions  used 
to  test  the  program  were  taken  from  Takeda  (1980,  p.  171) 
so  that  the  results  of  the  computer  program  could  be 
compared  with  experimental  results  from  Takeda  (1980)  and 
Chaturvedi  (1933).  The  properties  of  Scotchply  1002 
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composite  with  epoxy  matrix  and  E glass  fibers  used  in  the 
plate  are  as  follows: 

E1  = 40.0  GPa  (5.8  x 106  psi) 

E2  = 8.27  GPa  (1.2  x 106  psi) 

G 1 2 = 4.14  GPa  (0.6  x 10^  psi) 

9.1 

G23  = 3.03  GPa  (0.44  x 106  psi) 

v12  = 

V21  = 


Assuming  transverse  isotropy 


G13  - 212 


9.2 


The  density,  p,  with  units  of  mass  per  unit  volume 
determined  experimentally  was 


p = 1901.5  kg/m ^ 


9.3 


The  thickness,  h,  of  the  plate  was 


h = 0.00335  m 


9.4 


The  plate  was  clamped  on  the  edges  with  a length  and  width 
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of 


a = 0.14  m 


9.5 


between  the  inner  edges  of  the  supports.  This  corresponds 
to  the  value  given  by  Takeda  (1980),  but  because  the 
experimental  clamping  is  not  perfect,  the  effective 
experimental  dimension  may  be  slightly  larger. 

In  order  to  compare  analytical  solutions  for  isotropic 
plates  with  the  finite  element  results  and  experimental 
results  from  Takeda  (1980),  isotropic  material  properties 
which  approximate  the  anisotropic  properties  of  the  laminas 
were  chosen.  An  isotropic  Young’s  modulus  was  chosen  to  be 


which  yields 


E = 24.16  GPa 


9.7 


Poisson's  ratio  was  chosen  to  be 


v 


151 


9.8 


The  shear  modulus  was  calculated  by 


G 


E 


9.9 


2(1  + v) 
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which  yields 


= 10.469  GPa 


9.10 


Simply  Supported  Plate  with  a Concentrated 
Static  Load  at  the  Center 


The  finite  element  program  was  tested  by  solving  the 
problem  of  a simply  supported  isotropic  square  plate  with  a 
static  concentrated  load  at  the  center  as  shown  in  Figure 
9.1.  An  analytic  solution  for  classical  plate  theory  (CPT) 
for  a thin  plate  is  given  by  Timoshenko  and 
Woinowsky-Krieger  (1959,  p.  143) 


where  F is  the  concentrated  force,  a is  the  length  of  the 
side,  D is  the  bending  stiffness  of  the  plate  defined  in 
Equation  3.6,  and  wmax  is  the  displacement  of  the  plate  at 
the  center.  For  a square  plate 


Timoshenko  and  Woinowsky-Krieger  (1959,  p.  143)  give 


a 


9.12 


a = 0.01160 


9.13 


for  a square  plate,  which  is  correct  but  that  value  cannot 
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Figure  9.1. 


A Simply  Supported  Square  Plate  with  a 
Concentrated  Load  at  the  Center 
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be  obtained  from  Equation  9.7  as  written  in  Timoshenko  and 
Woinowsky-Krieger  (1959,  p.  143)  because  they  neglected  to 
say  that  the  series  must  be  summed  only  for  odd  terms. 

The  finite  element  program  was  used  to  calculate  the 
displacement  at  the  center  of  the  plate,  wmax,  which  was 
compared  with  displacements  calculated  using  Timoshenko's 
analytical  solution.  Equations  9.11  and  9.13.  The  program 
is  based  on  shear  deformable  plate  theory  which  included 
shear  deformation,  and  Timoshenko's  solution  is  based  on 
Classical  Plate  Theory  (CPT)  which  does  not  include  shear 
deformation.  A thin  plate  was  chosen  as  an  example  to 
minimize  the  shear  deformation  so  that  a valid  comparison 
between  the  results  of  the  two  methods  can  be  made.  The 
thickness  of  the  plate,  h,  and  the  length  of  the  side,  a, 
in  Equations  9.4  and  9.5  give  a ratio  of  length  to 
thickness,  a/h,  equal  42  which  results  in  small  shear 
deformation.  The  load,  F,  chosen  is  arbitrary  and  does  not 
affect  the  comparison  of  the  two  methods. 

In  order  to  non-diraensionalize  the  results,  Equation 
9.6  was  solved  for  <* 


a 


w D 
max 


Fa* 


9.14 


Values  of  a were  calculated  for  1,  4,  9,  and  25  elements 
and  the  results  are  shown  in  Table  9.1.  The  placement  of 
the  elements  on  the  plate  is  shown  in  Figure  3.1. 
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Table  9.1  Comparison  of  static  displacement  of  a simply 

supported  square  plate  with  a concentrated  load 
at  the  center. 


Values  of  a 


a = 


w D 
max 

Fa2 


number  of 

square  elements  SAP  IV 


Petersen 


1 

.01008 

.00829 

4 

.01095 

.01152 

9 

.01163 

25 

Ek— 1 

> 

.01163 

f 1 

T 0 O 

element  type 

< 

> o 

< 

5 

> < 

) O 6> 

Timoshenko  analytical  solution  a = .01160 
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The  values  of  a calculated  by  Petersen's  program 
converge  rapidly  towards  the  value  of  a from  Equation  9.11 
as  the  number  of  elements  increases.  For  a small  number  of 
elements  the  value  of  a calculated  by  the  computer  is  less 
than  the  analytical  value  of  a,  Equation  9.13,  because  a 
small  number  of  elements  cause  the  plate  to  be  overly 
stiff.  The  computed  value  of  a converges  to  a number 
slightly  larger  than  the  analytical  value  of  a because 
Timoshenko's  solution  is  based  on  CPT,  which  does  not 
include  shear  deformation  and  underestimates  the  deflection 
of  a plate.  The  agreement  was  judged  to  be  satisfactory. 

The  results  of  Petersen's  Program  were  also  compared 
with  results  from  SAP  IV,  a Structural  Analysis  program  for 
Static  and  Dynamic  Response  of  Linear  System  by  Bathe, 
Wilson,  and  Peterson  (1974).  The  same  problem,  a simply 
supported  isotropic  square  plate  with  a concentrated  load 
at  the  center,  was  used.  The  SAP  IV  program  is  based  on 
classical  plate  theory  (CPT)  and  uses  triangular 
elements.  Two  triangular  elements  are  used  to  obtain  each 
square  element.  The  results  are  shown  in  Table  9.1. 
Petersen's  program  converges  faster  than  SAP  IV  because 
each  square  element  in  Petersen's  program  contains  8 nodes 
whereas  each  square  element  in  SAP  IV  (two  triangles) 
contains  4 nodes,  and  the  SAP  IV  program  is  based  on  CPT 
which  under-predicts  deflections. 
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Cylindrical  Bending 

The  program  was  also  tested  by  solving  the  problem  of 
a plate  in  cylindrical  bending  as  shown  in  Figure  9.2. 
Since  the  plate  is  symmetric,  only  one  half  of  the  plate 
was  modeled.  Five  elements  were  used.  The  ends  were 
simply  supported.  Rotation  was  allowed  only  in  one 
direction,  which  means  the  plate  could  only  exhibit 
cylindrical  bending  plus  a small  shear  deformation.  A 
static  force  perpendicular  to  the  plate  was  applied  at  the 
center.  Cylindrical  bending  of  a plate  was  chosen  because 
it  is  similar  to  beam  deflection,  for  which  there  are 
analytical  solutions. 

Cylindrical  bending  of  an  isotropic  plate  is  governed 
by 


D w 


XX 


M 


where  M is  the  moment  and 


9.15 


Eh 

2(1  - v2) 


9.16 


If  v is  zero 


D = E I/b 


9.17 
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Figure  9.2.  Plate  in  Cylindrical  Bending 
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where  b is  the  width  of  the  plate  and  the  moment  of 
inertia,  I,  is 


I 


9.18 


and  Equation  9.15  becomes  the  same  as  the  differential 
equation  of  a beam  of  unit  width.  The  program  was  run  with 
v equal  zero,  and  the  results  were  compared  with  beam 
theory. 

The  maximum  normal  stress,  CTmax,  in  an  isotropic  plate 
in  cylindrical  bending  occurs  at  the  top  and  bottom 
surfaces  and  can  be  calculated  for  v = 0 using  the  equation 
for  a beam 


a = MJL 
max  I 2 


9.19 


The  finite  element  program  calculates  the  normal  stress 
with  accuracy  of  at  least  5 digits  compared  with  results 
from  Equation  9.19.  The  average  shear  stress,  Tav»  in  a 
narrow  isotropic  plate  in  cylindrical  bending  shown  in 
Figure  9.1  can  be  calculated  for  points  not  near  the 
concentrated  load  using 


T 


_ F 

av  _ bh 


9.20 


The  finite  element  program  calculates  the  average  shear 
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stress  with  accuracy  of  at  least  5 digits  compared  with 
results  from  Equation  9.20. 

The  deflection  of  the  midpoint  of  a simply  supported 
beam  with  load  at  the  center  according  to  Thompson  (1965, 
p.  25)  is 


w 


FL' 


48EI 


9.21 


where  L is  the  length  of  the  beam.  For  a thin  beam  with 
L/h  ratio  of  42  the  difference  between  the  deflections 
computed  by  Petersen's  program  and  the  deflections  computed 
by  Equation  9.21  was  less  than  0.2  percent. 


Period  of  Free  Vibration 

The  program  was  tested  by  computing  the  period  of  free 
vibration  and  comparing  this  with  an  analytical  solution 
and  with  experimental  results.  The  period  of  free 
vibration  of  an  isotropic  plate  is  given  in  Timoshenko, 
Young,  and  Weaver  (1974,  p.  494).  For  a simply  supported, 
isotropic,  square  plate  the  period  of  the  lowest  mode  of 
vibration  is 

2 

P1  = f-  7 (D/ph)  9.22 

The  period  of  the  lowest  mode  was  calculated  using 
isotropic  material  properties,  Equation  9.6,  density, 
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Equation  9.3»  and  plate  dimensions,  Equations  9.4  and 
9*5.  The  period  is  shown  in  Table  9*2.  Petersen’s  program 
was  tested  by  calculating  the  period  for  the  same  mode  as 
given  by  the  analytical  solution.  The  same  isotropic 
material  properties  and  plate  dimensions  were  used.  The 
plate  was  given  an  initial  velocity.  The  deflection  of  the 
plate  was  calculated  for  each  time  step.  The  period  of 
first  mode  of  variation  was  identified  as  the  time  when  the 
plate  returned  to  its  initial  position. 

The  lowest  period  of  free  vibration  of  a simply 
supported  square  isotropic  plate  calculated  by  Petersen's 
program  for  1,  9,  and  25  elements  is  shown  in  Table  9.2. 

The  period  calculated  by  Petersen's  program  converges 
quickly  toward  the  period  calculated  by  Equation  9.22  as 
the  number  of  elements  is  increased.  The  period  calculated 
analytically  is  a little  greater  than  the  period  calculated 
by  Petersen's  program  because  the  analytical  solution  is 
based  on  CPT  which  neglects  rotary  inertia  and  shear 
deformation  and  overestimates  the  period.  The  period  for 
the  same  material  properties  and  plate  dimension  calculated 
using  SAPIV  program  based  on  CPT  is  shown  in  Table  9.2  for 
comparison . 
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Table  9.2  Lowest  period  of  vibration  of  an  isotropic 
simply  supported  square  plate. 


microseconds 


number  of 
square  elements 

SAP  IV 

Petersen 

1 

1450 

4 

1760 

9 

1750 

25 

1765 

analytical  solution  1790 


CHAPTER  TEN 

DISPLACEMENTS  AND  STRESSES  IN  A PLATE 
IMPACTED  BY  A PROJECTILE 

Introduction 

The  finite  element  program  developed  for  this  study 
was  used  to  solve  the  problem  of  elastic  waves  and 
vibrations  in  a square  plate  subjected  to  loadings 
simulating  impact  by  a projectile  at  the  center.  The  plate 
and  projectile  are  shown  in  Figure  10.1.  The  x and  y 
coordinate  axes  are  in  the  plane  of  the  plate.'  The 
projectile  moves  in  the  positive  z direction,  which 
produces  a positive  initial  deflection  at  the  center  of  the 
plate.  The  stresses  were  calculated  at  the  Gauss  points, 
which  are  identified  in  Figure  10.1  with  x's  near  the  x or 
y axes.  There  are  two  Gauss  points  near  the  x axis  for 
each  element  bordering  the  x axis. 

The  plate  materials  may  be  isotropic,  orthotropic,  or 
layered  orthotropic.  The  impact  can  be  either  elastic  or 
inelastic.  The  boundary  conditions  can  be  clamped  or 
simply  supported.  The  number  of  elements  and  the  time  step 
can  be  chosen  to  produce  the  desired  accuracy. 

When  a projectile  impacts  a plate,  waves  are 
produced.  In  this  chapter  the  deflections,  shear  stresses, 
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Figure  10.1.  Projectile  and  Plate  with  Gauss  Points 
Indicated 
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and  normal  stresses  are  displayed  graphically  so  that  the 
waves  can  be  observed.  The  waves  are  important  because 
experiments  conducted  by  Takeda,  Sierakowski,  Ross,  and 
Malvern  (1982)  indicate  that  flexural  wave  propagation  is 
the  principal  response  mode  of  moderately  thin  and  thick 
plates  to  central  impact,  and  it  leads  to  delamination- 
crack  propagation.  The  analysis  here  is  concerned  only 
with  the  elastic  wave  propagation,  assuming  that 
delamination  has  not  yet  occurred.  The  resulting  stress 
distributions  indicate  where  and  when  delamination  might 
begin. 


Isotropic  Plate  with  Inelastic  Impact 
The  results  for  an  isotropic  plate  will  be  presented 
first.  The  plate  was  assumed  to  have  the  isotropic 
material  properties  as  given  in  Equations  9.7,  9.8,  and 
9-10,  density  as  given  in  Equation  9.3»  and  dimensions  as 
given  in  Equations  9.4  and  9.5.  The  irapactor  has  mass 
0.014175  kg  and  velocity  19.8  m/second.  The  plate  is 
clamped  with  a side  length  of  0.14  m.  These  dimensions  and 
impactor  mass  are  the  same  as  those  used  in  the  experiments 
of  Takeda  (1980),  and  the  isotropic  properties  assumed  are 
averages  of  the  orthotropic  properties  he  assumed.  In  this 
first  analysis  the  impact  is  assumed  to  be  inelastic,  with 
the  projectile  mass  being  added  to  the  mass  of  the  impacted 
node  for  subsequent  analysis  of  plate  motion. 
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The  deflection  of  a plate  impacted  by  a projectile  is 
shown  in  Figure  10.2.  A pictorial  drawing  of  the 
deflection  of  the  plate  is  shown  in  Figure  10.2a.  The 
abscissa  is  the  x axis  in  the  plane  of  the  undeformed 
plate,  and  the  ordinate  is  the  deflection,  w,  of  the 
plate.  The  plate  shown  in  Figure  10.2a  has  been  cut  along 
the  x axis  to  show  a cross-sectional  view  of  the  deflected 
plate.  The  four  views  shown  are  for  successive  times  as 
the  projectile  deflects  the  plate.  The  deflection  of  the 
plate  was  calculated  by  the  finite  element  computer  program 
and  displayed  graphically  in  Figure  10.2b.  Since  the  plate 
is  symmetric  with  respect  to  the  x and  y axes,  only  one 
quarter  of  the  plate  was  analyzed  by  finite  elements.  The 
plate  deflection  is  shown  for  five  successive  times. 

An  enlarged  view  of  the  calculated  deflection  of  the 
isotropic  plate  hit  by  a projectile  is  shown  in  Figure 
10.3*  The  deflection  is  measured  from  the  static 
equilibrium  position.  The  deflection  along  the  x axis  is 
shown.  The  x axis  gives  the  distance  from  the  center  of 
the  plate  where  impact  occurs.  The  center  of  the  plate  is 
at  x equal  zero,  and  the  edge  of  the  plate  is  at  x equal  70 
mm.  The  deflection  is  shown  for  five  different  times:  10, 

20,  30,  40,  and  50  microseconds,  which  are  identified  by 
the  legend. 

The  plate  is  initially  at  rest  with  no  deflection  at 
time  equal  0.  At  time  equal  10  microseconds  the  plate  has 
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Figure  10.2.  Cross  Section  of  a Plate  Impacted  by  a 
Projectile 

a)  Pictorial  Deflection,  b)  Computed 
Deflection 
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deflected  near  the  center,  but  remains  stationary  near  the 
edges.  It  takes  a small  but  finite  time  for  the 
disturbance  at  the  center  to  propagate  to  the  edges.  The 
deflection  of  the  center  of  the  plate  increases  as  time 
increases  from  0 to  50  microseconds,  and  the  flexural  wave 
propagates  toward  the  boundary. 

The  effects  of  boundary  conditions  can  be  observed  in 
Figure  10.3.  The  slope  is  zero  at  x equal  zero  because  the 
x and  y axes  are  lines  of  symmetry.  The  slope  and 
deflection  are  zero  at  x equal  70  ram  because  the  edges  of 
the  plate  are  clamped. 

The  deflection  is  positive  near  the  center  of  the 
plate  but  negative  at  some  other  places.  The  negative 
deflection  is  caused  by  the  bending  stiffness  and  the  mass 
of  the  plate.  Similar  negative  deflection  can  be  seen  in 
the  response  of  an  infinite  plate  to  a concentrated  load  at 
the  origin  presented  by  Meirovitch  (1967,  p.  381).  At 
longer  times,  after  reflections  of  the  flexural  wave  from 
the  boundary,  negative  deflections  would  also  appear  at  the 
center  as  the  plate  goes  into  vibration. 

The  normal  stress,  ax,  on  the  surface  of  the  plate  is 
shown  in  Figure  10.4.  The  stress  is  calculated  on  the 
surface  where  impact  occurs  (z  equal  -h/2)  which  gives  the 
maximum  or  minimum  normal  stress.  The  stress  is  calculated 
at  the  Gauss  points  near  the  x axis,  which  are  shown  in 
Figure  10.1.  The  stress  is  not  calculated  exactly  on  the 
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x axis  because  stress  calculations  are  more  accurate  at  the 
Gauss  points  according  to  Zienkiewicz  (1977,  p.  281). 

The  largest  magnitude  normal  stress  is  a compressive 
stress  near  the  point  of  impact.  Stress  near  the  point  of 
impact  cannot  be  calculated  accurately  by  this  finite 
element  program.  The  maximum  tensile  normal  stress  can  be 
observed  in  Figure  10.4  at  10  mm  from  the  center  of  the 
plate  for  time  equal  10  microseconds.  The  position  of  the 
maximum  tensile  normal  stress  moves  from  the  center  towards 
the  edge  of  the  plate.  The  motion  of  the  point  of  maximum 
tensile  stress  is  caused  by  a large  flexural  wave  generated 
by  the  initial  impact.  The  flexural  wave  begins  at  the 
center  of  the  plate  and  moves  toward  the  edge. 

The  shear  stress,  Txz,  averaged  through  the  thickness 
of  the  plate,  is  shown  in  Figure  10.5.  The  shear  stress 
was  calculated  at  the  Gauss  points  near  the  x axis  where 
the  normal  stress  was  calculated.  The  largest  magnitude 
shear  stress  occurs  near  the  impact  point,  but  this  stress 
cannot  be  calculated  accurately  by  the  finite  element 
program.  The  position  of  maximum  positive  shear  stress 
moves  from  the  center  of  the  plate  towards  the  edge  with 
increasing  time.  The  motion  of  the  maximum  positive  shear 
stress  is  caused  by  a flexural  wave  moving  from  the  point 
of  impact  toward  the  edge. 
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Figure  10.3*  Deflection  of  an  Isotropic  Plate  for  Time 
from  10  to  50  microseconds 
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Figure  10.4.  Normal  Stress,  crx,  in  an  Isotropic  Plate  for 
Time  from  10  to  50  microseconds 
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Figure  10.5.  Average  Shear  Stress,  Txz»  in  an  Isotropic 
Plate  for  Time  from  10  to  50  microseconds 
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Fifteen-Layer  Plate  with  Inelastic  Impact 

The  finite  element  program  will  calculate  deflections 
and  stresses  for  layered  orthotropic  materials  as  well  as 
for  isotropic  materials.  A particular  15-layered  material 
was  chosen  so  that  results  of  the  finite  element  analysis 
could  be  compared  with  the  deflections  determined 
experimentally  by  Chaturvedi  (1983).  The  plate  was  made 
of  glass-epoxy  laminas  with  material  properties  as 
given  in  Equations  9.1,  9.2,  and  9.3  and  a layup 
[0°/90o/0°  . . . 1-J5*  The  plate  dimensions  are  given  in 
Equations  9.4  and  9.5.  The  plate  was  clamped,  but  it  is 
impossible  to  produce  a truly  clamped  boundary  condition  in 
an  experiment.  The  plate  was  impacted  with  a projectile 
with  mass  0.014175  kg  and  velocity  19.8  m/s.  The  impact 
was  assumed  to  be  inelastic,  the  same  as  for  the  isotropic 
plate . 

The  curves  showing  the  deflection  of  the  center  of  the 
plate  as  a function  of  time  determined  experimentally  and 
by  the  finite  element  program  are  shown  in  Figure  10.6. 

The  curves  showing  the  deflection  of  the  center  of  the 
plate  calculated  by  the  finite  element  program  are  for  an 
inelastic  collision.  One  curve  is  for  a simply  supported 
plate,  and  the  other  is  for  a clamped  plate.  The  clamped 
plate  has  a shorter  period  than  the  simply  supported 
plate.  The  deflection  calculated  with  clamped  boundary 
conditions  agrees  with  the  deflection  determined 
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Experiment  by  Chaturvedi  (1983) 
Finite  Element  Analysis  with  Simply 
Supported  Boundary  Conditions 
Finite  Element  Analysis  with  Clamped 
Boundary  Conditions 


Figure  10.6.  Deflection  of  the  Center  of  a 15-Layer  Plate 
Determined  Experimentally  and  by  Finite 
Element  Analysis 
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experimentally  for  time  less  than  400  microseconds.  As  the 
deflection  of  the  plate  approaches  the  maximum  at  500 
microseconds,  the  experimental  and  computed  deflections  do 
not  agree  because  the  experimental  boundary  condition  is 
not  truly  clamped. 

The  metal  frame  cannot  hold  the  impacted  plate  tightly 
enough  to  produce  a truly  clamped  boundary,  where  the  plate 
deflection  and  rotation  are  exactly  zero  at  the  boundary. 
The  metal  frame  allows  the  plate  to  rotate  a small  amount, 
which  is  neither  a truly  clamped  nor  a simply  supported 
boundary,  but  somewhere  between  the  two  conditions. 

The  deflection,  normal  stress,  ax,  on  the  surface  (z 
equal  -h/2),  and  average  shear  stress,  txz,  in  a 15-layered 
plate  near  the  x axis  for  time  from  10  to  50  microseconds 
are  shown  in  Figures  10.7,  10.8,  and  10.9.  The  shear 
stress  is  an  average  shear  stress,  which  is  found  by  taking 
the  arithmetic  mean  of  the  shear  stresses  in  the  various 
layers.  The  deflections  and  shear  stresses  are  similar  in 
magnitude  to  the  deflections  and  shear  stresses  shown  in 
Figures  10.3  and  10.5  for  an  isotropic  plate  of  similar 
materials,  and  the  same  size,  boundary  conditions,  and 
loading.  The  normal  stress,  ax,  is  larger  in  the  15-layer 
plate  shown  in  Figure  10.3  than  in  the  isotropic  plate 
shown  in  Figure  10.4  for  the  reason  which  follows.  The 
normal  stress  is  calculated  on  the  outside  surface  of  the 
plate.  For  the  15-layer  plate,  the  modulus  of 
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Deflection  of  a 15-Layer  Plate  Along  the  x 
Axis  for  Time  from  10  to  50  microseconds 


Figure  10.7. 
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Figure  10.8.  Normal  Stress,  a„,  in  a 15-Layer  Plate  Along 
the  x Axis  for  Time  from  10  to  50 
microseconds 
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Figure  10.9.  Average  Shear  Stress,  t , in  a 15-Layer 

Plate  Along  the  x Axis  for  Time  from  10  to  50 
microseconds 
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elasticity  in  the  outside  layer  in  the  x direction  is 
larger  than  the  modulus  of  elasticity  in  the  isotropic 
plate. 

The  deflection,  normal  stress,  o , and  average  shear 
stress,  TyZ*  near  the  y axis  in  a 15-layer  plate  are  shown 
in  Figures  10.10,  10.11,  and  10.12  for  time  from  10  to  50 
microseconds.  The  bending  stiffness  of  the  plate  in  the  x 
direction,  D^,  is  greater  than  the  bending  stiffness  in 
the  y direction,  D22’  for  a plate  with  layup 
[0°/90°/0°  . . . ]15  and  E1  greater  than  E£.  A higher 
bending  stiffness  in  the  x direction  than  in  the  y 
direction  causes  the  flexural  waves  to  move  faster  in  the  x 
direction  than  the  y direction.  This  can  be  observed  by 
comparing  Figures  10.7,  10.8,  and  10.9  with  Figures  10.10, 
10.11,  and  10.12.  The  points  of  maximum  normal  stress  and 
shear  stress  move  faster  in  the  x direction  than  in  the  y 
direction  because  the  flexural  waves  move  faster  in  the  x 
direction  than  in  the  y direction.  The  normal  stress,  cjx, 
is  greater  near  the  x axis  than  the  normal  stress,  Oy,  near 
the  y axis  because  the  modulus  of  elasticity  in  the  surface 
layer  is  greater  in  the  x direction  than  in  the  y 
direction . 

The  deflection,  normal  stress,  oy,  and  average  shear 
stress,  Txz,  along  the  x axis  in  a 15-layer  plate  are  shown 
in  Figures  10.13,  10.14,  and  10.15  for  time  from  100  to  500 
microseconds.  The  large  flexural  wave  is  approaching  the 
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Figure  10.10.  Deflection  of  a 15-Layer  Plate  Along  the  y 

Axis  for  Time  from  10  to  50  microseconds 
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Figure  10.11.  Normal  Stress,  a , in  a 15-Layer  Plate  Along 

the  y Axis  for  Time  from  10  to  50 
microseconds 
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Figure  10.12.  Average  Shear  Stress,  t„z,  in  a 15-Layer 

Plate  Along  the  y Axis  ror  Time  from  10  to 
50  microseconds 
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Figure  10.13*  Deflection  of  a 15-Layer  Plate  Along  the  x 

Axis  for  Time  from  100  to  500  microseconds 
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Figure  10.14.  Normal  Stress,  ax>  in  a 15-Layer  Plate  Along 

the  x Axis  for  Time  from  100  to  500 
microseconds 


127 


LEGEND:  TIME  e-e-B  100  a- a-  a 20D 
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Figure  10.15*  Average  Shear  Stress,  Txz»  a “15-Layer 

Plate  Along  the  x Axis  for  Time  from  100  to 
500  microseconds 


128 


boundary  at  time  equal  100  microseconds,  producing  large 
stresses  at  the  boundary.  The  deflection  increases  after 
time  equals  100  microseconds  but  the  stresses  do  not  in 
general  increase.  The  deflection  of  the  center  of  the 
plate  reaches  a maximum  between  400  and  500  microseconds 
(see  Figure  10.6),  but  the  stresses  are  not  maximum  at  the 
same  time. 

The  boundary  conditions  were  changed  from  clamped  to 
simply  supported  to  investigate  the  effect  of  boundary 
conditions  on  the  plate.  The  deflection,  normal  stress, 
ax,  and  average  shear  stress,  xxz*  f°r  simply  supported 
boundary  conditions  at  time  from  10  to  50  microseconds  were 
almost  the  same  as  the  stresses  shown  in  Figures  10.7, 

10.8,  and  10.9  for  clamped  boundary  conditions.  The  waves 
reflected  from  the  boundary  before  50  microseconds  are  not 
significant . 

The  deflections  and  stresses  in  a 15-layer  simply 
supported  plate  are  shown  in  Figures  10.16,  10.17,  and 
10.18  for  time  from  100  to  500  microseconds.  The 
differences  between  the  deflections  and  stresses  for  simply 
supported  boundary  conditions  and  clamped  boundary 
conditions,  Figures  10.13,  10.14,  and  10.15,  are  small  at 
time  equal  100  microseconds.  The  differences  increase  with 
time  as  large  flexural  waves  are  reflected  from  the 


boundaries . 


129 


LEGEND:  TIME  □ □ Q 200  a A a 20D 

# #•  » 300  h — i — >-  ^00 

z-  z-g  500 


Figure  10.16.  Deflection  of  a 15-Layer  Plate  Along  the  x 

Axis  for  Time  from  100  to  500  microseconds 
for  Simply  Supported  Boundary  Conditions 
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Figure  10.17.  Normal  Stress,  o in  a 15-Layer  Plate  Along 

the  x Axis  for  Time  from  100  to  500 
microseconds  for  Simply  Supported  Boundary 
Conditions 
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Figure  1 
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Average  Shear  Stress,  T„z,  in  a 15-Layer 
Plate  Along  the  x Axis  for  Time  from  100  to 
500  microseconds  for  Simply  Supported 
Boundary  Conditions 
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The  deflection  of  the  center  of  the  plate  at  time 
equal  400  microseconds  is  greater  for  the  plate  with 
clamped  boundary  conditions  than  for  the  plate  with  simply 
supported  boundary  conditions  as  can  be  seen  in  Figure 
10.6.  The  relatively  small  center  deflection  can  be 
explained  by  observing  the  deflection  of  a simply  supported 
plate  at  time  equal  400  microseconds  shown  in  Figure 
10.16.  The  deflection  of  the  plate  at  the  center  is  less 
than  the  deflection  near  the  center. 

The  plate  thickness,  the  projectile  velocity,  and  the 
projectile  mass  were  altered  to  determine  the  effects  these 
factors  have  on  the  deflections  and  stresses.  The  velocity 
of  the  impactor  was  doubled;  it  was  increased  from  19.8 
m/second  to  38.6  m/second  and  the  results  are  shown  in 
Figures  10.19,  10.20,  and  10.21.  The  only  change  is  that 
the  magnitudes  of  the  deflections  and  stresses  are  doubled 
as  compared  with  the  results  for  a velocity  of  19.8 
m/second  shown  in  Figures  10.7,  10.8,  and  10.9.  The 
magnitudes  of  the  deflections  and  stresses  vary  linearly 
with  the  magnitude  of  the  projectile  velocity,  provided  the 
response  remains  linearly  elastic,  as  was  assumed  for  the 
finite  element  analysis. 

The  computed  deflection,  normal  stress,  ax,  and  average 
shear  stress,  Txz»  along  the  x axis  for  a 15-layer  plate  im- 
pacted by  a projectile  with  mass  0.14175  kg  at  time  from  10  to 
50  microseconds  are  shown  in  Figures  10.22,  10.23,  and  10.24. 
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Figure  10.19.  Deflection  of  a 15-Layer  Plate  Along  the  x 

Axis  for  Time  from  10  to  50  microseconds  for 
Impactor  Velocity  of  58.6  m/s 
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Figure  10.20.  Normal  Stress,  cr^,  in  a 15-Layer  Plate  Along 

the  x Axis  for  Time  from  10  to  50  micro- 
seconds for  Impactor  Velocity  of  38.6  m/s 
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Figure  10.21.  Average  Shear  Stress,  txz,  in  a 15-Layer 

Plate  Along  the  x Axis  for  Time  from  10  to 
50  microseconds  for  Impactor  Velocity  of 
38.6  m/s 
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Figure  10.22.  Deflection  of  a 15-Layer  Plate  Along  the  x 

Axis  for  Time  from  10  to  50  microseconds  for 
Projectile  with  Mass  0.14175  kg 
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Figure  10.23.  Normal  Stress,  a in  a 15-Layer  Plate  Along 

the  x Axis  for  Time  from  10  to  50  micro- 
seconds for  Projectile  with  Mass  0.14175  kg 
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Figure  10.24.  Average  Shear  Stress,  t , in  a 15-Layer 

Plate  Along  the  x Axis  for  Time  from  10  to 
50  microseconds  for  Projectile  with  Mass 
0.14175  kg 
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The  mass  of  the  projectile  is  10  times  larger  than  the  mass 
of  the  projectile  used  in  the  calculations  for  Figures 
10.7,  10.8,  and  10.9.  The  deflections  and  stresses  are 
increased  less  than  15  percent  when  the  mass  of  the 
projectile  is  increased  by  a factor  10.  This  means  that 
for  time  less  than  50  microseconds  an  increase  in  the  mass 
of  the  projectile  has  little  effect  on  the  deflections  and 
stresses  in  the  plate. 

The  computed  deflection,  normal  stress,  ax,  and  shear 
stress,  fx2,  along  the  x axis  of  a 15-layer  plate  with 
total  thickness  0.0067  m for  time  from  10  to  50  micro- 
seconds are  shown  in  Figures  10.25,  10.26,  and  10.27.  The 
projectile  had  a mass  of  0.014175  kg  and"  velocity  of  19.8 
m/second.  The  15-layer  plate  is  twice  as  thick  as  the 
plate  used  for  the  calculated  deflections  and  stresses  in 
Figures  10.7,  10.8,  and  10.9.  When  the  thickness  of  the 
plate  is  doubled,  the  deflection  and  normal  stress,  crx,  are 
reduced  a small  amount  but  the  shear  stresses  show  some 
increases  and  some  decreases.  The  flexural  wave  velocity 
is  greater  for  the  thicker  plate.  The  points  of  maximum 
normal  and  shear  stress  move  toward  the  edges  faster  in  the 
thicker  plate  than  in  the  thinner  plate. 

Three-Layer  Plate  with  Inelastic  Impact 

In  the  previous  section,  the  deflections  and  stresses 
in  a plate  with  15-layers  and  a projectile  velocity  of 
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Figure  10.25.  Deflection  of  a 15-Layer  Plate  Along  the  x 

Axis  for  Time  from  10  to  50  microseconds  for 
Plate  Thickness  0.0067  m 
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Figure  10.26.  Normal  Stress,  in  a 15-Layer  Plate  Along 

the  x Axis  for  Time  from  10  to  50 
microseconds  for  Plate  Thickness  0.0067  m 
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Figure  10.27.  Average  Shear  Stress,  x„z,  in  a 15-Layer 

Plate  Along  the  x Axis  for  Time  from  10  to 
50  microseconds  for  Plate  Thickness  0.0067  ni 
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19.8  m/second  were  presented.  In  this  section,  the 
deflections  and  stresses  in  a plate  with  3-layers  and  a 
projectile  velocity  of  37.3  m/second  will  be  presented. 

The  same  inelastic  impact  law  was  used.  The  plate  had  the 
same  material  properties  as  given  in  Equations  9.1,  9.2, 
and  9.3  and  plate  dimensions  as  given  in  Equations  9.4  and 
9.5  as  the  previous  plate.  The  layup  was  [0° ^ / 90° ^ / 0° ^ ] . 
This  plate  was  used  for  experimental  work  by  Takeda  (1980), 

Takeda,  Sierakowski,  and  Malvern  (1981a;  1981b),  and 

* 

Takeda,  Sierakowski,  Ross,  and  Malvern  (1982). 

The  deflection,  normal  stress,  °x,  and  the  average 
shear  stress,  Txz,  along  the  x axis  in  the  3-layer  plate 
are  shown  in  Figures  10.28,  10.29,  and  10.30  for  time  from 
10  to  50  microseconds  after  impact  and  in  Figures  10.31, 
10.32,  and  10.33  for  time  from  50  to  250  microseconds.  The 
deflection,  normal  stress,  ay,  and  the  average  shear 
stress,  T,,_,  along  the  y axis  in  the  3-layer  plate  for  the 
same  impact  are  shown  in  Figures  10.34,  10.35,  and  10.36 
for  time  from  10  to  50  microseconds  and  in  Figures  10.37, 
10.38,  and  10.39  for  time  from  50  to  250  microseconds. 

The  bending  stiffness  is  greater  than  the  bending 

stiffness  D22  for  a 3-layer  plate  with  layup  [0°/90°/0°], 
and  E-j  is  greater  than  E2.  This  difference  in  bending 
stiffness  causes  the  flexural  wave  to  move  faster  in  the  x 
direction  than  in  the  y direction.  This  difference  in 
flexural  wave  velocity  can  be  observed  in  the  deflections 
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Figure  10.28.  Deflection  of  a 3-Layer  Plate  Along  the  x 

Axis  for  Time  from  10  to  50  microseconds 
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Figure  10.29.  Normal  Stress,  crx,  in  a 3-Layer  Plate  Along 

the  x Axis  for  Time  from  10  to  50 
microseconds 
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Figure  10.31.  Deflection  of  a 3-Layer  Plate  Along  the  x 

Axis  for  Time  from  50  to  250  microseconds 
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Figure  10.32.  Normal  Stress,  ox,  in  a 3-Layer  Plate  Along 

the  x Axis  for  Time  from  50  to  250 
microseconds 


149 


LEGEND:  TIME  BOB  50  * ft  & 10D 

* » » 150  i — i — i-  200 

z--z -z  250 


Figure  10.33.  Average  Shear  Stress,  t in  a 3-Layer 

Plate  Along  the  x Axis  ror  Time  From  50  to 
250  microseconds 
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Figure  10.34.  Deflection  of  a 3-Layer  Plate  Along  the  y 

Axis  for  Time  from  10  to  50  microseconds 
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Figure  10.35.  Normal  Stress,  ay,  in  a 3-Layer  Plate  Along 

the  y Axis  for  Time  from  10  to  50 
microseconds 
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Figure  10.56.  Average  Shear  Stress,  t , in  a 5-Layer 

Plate  Along  the  y Axis  ror  Time  from  10  to 
50  microseconds 
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Figure  10.37*  Deflection  of  a 3-Layer  Plate  Along  the  y 

Axis  for  Time  from  50  to  250  microseconds 
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Figure  10.38.  Normal  Stress,  ay,  in  a 3-Layer  Plate  Along 
the  y Axis  for  Time  from  50  to  250 
microseconds 
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Figure  10.39. 


250  microseconds 
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and  stresses.  For  example,  observe  the  plate  deflection  at 
time  equal  50  microseconds  shown  in  Figures  10.28  and 
10.34.  The  distances  from  the  origin  to  the  nearest  point 
where  the  plate  intersects  the  x or  y axis  is  30  mm  in  the 
x-direction  (Figure  10.28)  and  22  mm  in  the  y-direction 
(Figure  10.34).  The  3-layer  plate  exhibits  greater 
anisotropy  than  the  15-layer  plate  of  the  same  thickness. 

Elastic  Impact 

The  deflections  and  stresses  in  the  previous  section 
were  calculated  with  the  assumption  of  an  inelastic 
impact.  In  this  section  an  elastic  impact  was  simulated. 
The  elastic  impact  allows  the  projectile  to  indent  the 
plate  and  rebound.  The  projectile  is  allowed  to  lose 
contact  with  the  plate,  which  was  not  possible  for  the 
inelastic  collision.  Therefore,  the  elastic  impact  gives 
more  accurate  displacement  and  stress  after  about  350 
microseconds  when  the  projectile  loses  contact  with  the 
plate . 

The  force,  F,  between  the  projectile  and  the  plate 
during  impact  was  assumed  to  be  governed  by  the  impact  law 

F = H (r  - w)  p 7.14 

where  H is  a material  constant,  r is  the  displacement  of 
the  projectile,  w is  the  displacement  of  the  plate,  and  p 
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is  the  exponent,  as  discussed  in  Chapter  Seven.  For  p 
equal  1,  the  relation  between  force  and  displacement  is 
linear.  For  p equal  1.5,  Equation  7.3  is  the  Hertzian 
impact  law.  Both  the  linear  and  the  Hertzian  impact  law 
were  tried.  An  attempt  was  made  to  choose  a value  of  H 
which  would  cause  the  computed  deflections  to  match  the 
experimentally  determined  deflections  shown  in  Figure  10.6. 

A linear  impact  law  with  p equal  1 was  tried  and 
abandoned.  Different  values  of  H were  tried,  but  no  value 
could  be  found  to  produce  computed  deflections  which 
matched  the  experimentally  determined  deflections.  Also 
the  computed  deflections  were  sensitive  to  the  value  of  H; 
a small  change  in  H produced  a large  change  in  the 
deflection . 

Fifteen-Layer  Plate  with  Hertzian  Impact 

A Hertzian  impact  law  (p  equal  1.5)  was  tried  and  used 

successfully.  Various  values  of  H were  tried.  The 

computed  and  experimentally  determined  deflections  of  the 

center  of  a 15-layer  plate  as  a function  of  time  are  shown 

in  Figure  10.40.  The  15-layer  plate  had  the  same  material 

properties,  thickness,  as  the  previous  15-layer  plate  and 

the  projectile  velocity  was  the  same,  19-8  m/seconds.  The 

length  of  the  side  was  chosen  to  be  0.16  m,  and  value  of  H 

6 IS 

was  chosen  to  be  100  x 10  N/m  . 
explained  later. 


These  choices  will  be 
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Experiment  by  Chaturvedi  (1983) 

Finite  element  analysis 

Projectile  in  contact  with  plate 

Figure  10. 40.  Deflection  of  the  Center  of  a 15-Layer  Plate 

Computed  with  H Equal  100  x 10^  and 

Side  Length  0.16  ra,  for  Hertzian  Impact 
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The  agreement  between  the  deflections  determined 
experimentally  and  deflections  computed  using  a Hertzian 
law  shown  in  Figure  10.40  is  almost  as  good  as  the  excel- 
lent agreement  between  the  deflections  determined  experi- 
mentally and  the  deflections  computed  using  an  inelastic 
collision  shown  in  Figure  10.6  for  time  less  than  400 
microseconds.  For  time  greater  than  400  microseconds  the 
Hertzian  impact  law  gives  more  accurate  deflections  because 
the  Hertzian  law  allows  the  projectile  to  leave  the  plate. 

The  computed  deflection  of  the  center  of  a 15-layer 
plate  are  shown  in  Figure  10.41  for  H equal  50  x 10^,  100  x 
10^,  and  150  x 10^  N/m"'*^.  Note  that  the  positive  axis  is 
down  in  Figure  10.40  and  up  in  Figure  10.41.  Figure  10.41 
was  produced  using  fewer  elements  and  larger  time  steps 
than  Figure  10.40  and  the  results  may  show  differences. 

The  deflections  are  the  same  for  the  three  values  of  H 
shown,  which  means  the  computed  deflections  are  not 
sensitive  to  the  value  of  H chosen.  The  Hertzian  impact 
law  is  a good  choice  because  a precise  value  of  H is  not 
needed.  The  value  of  H has  been  determined  experimentally 
by  Sun,  Sankar,  and  Tan  (1981)  for  static  loading,  but  the 
value  of  H for  dynamic  loading  and  unloading  is  unknown. 

The  forces  between  projectile  and  the  plate  for  H 
equal  50  x 10^  N/m1*^  and  H equal  150  x 10  N/m1*^  are  shown 
in  Figure  10.42,  and  for  H equal  100  x 10^  N/m^'^  are  shown 
in  Figure  10.43.  Figure  10.43  was  computed  using  more 
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Figure  10.42.  Force  Between  the  Projectile  and  a 15-Layer  Plate  for  H Equal 
50  x 10^  N/m1*^  and  for  H Equal  150  x 10^  N/m^*^ 
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elements  and  smaller  time  steps  than  Figure  10.42,  and  the 
results  may  show  small  differences.  Note  that  the  smaller 
value  of  H,  which  models  a softer  material,  produces  a 
smoother  curve  of  force  versus  time.  Since  the  maximum 
deflections  of  the  plate  in  Figure  10.41  are  almost  the 
same,  the  total  area  under  the  two  curves  in  Figure  10.42 
should  be  the  same.  Note  that  the  magnitude  of  the  force 
computed  with  the  smaller  value  of  H lags  slightly  behind 
the  magnitude  of  the  force  computed  with  the  larger  value 
of  H.  The  force  computed  with  H equal  150  x 10^  N/m1*^ 
shows  some  numerical  instability,  which  can  be  also 
observed  in  the  deflections  shown  in  Figure  10.41c.  For 
values  of  H larger  than  150  x 10^  N/m^*^  numerical 
instabilities  make  the  results  unacceptable. 

When  the  force  remains  at  zero,  the  projectile  has 
lost  contact  with  the  plate.  The  first  impact  occurs  from 
0 to  approximately  350  microseconds,  at  which  time,  the 
projectile  separates  from  the  plate.  The  projectile 
recontacts  the  plate  at  approximately  700  microseconds. 

The  projectile  recontacts  the  plate  a third  time  at 
approximately  950  microseconds. 

The  multiple  contacts  can  be  observed  in  the  experi- 
mentally determined  curve  in  Figure  10.40.  The  center  of 
the  plate  accelerates  downward  at  300  microseconds  causing 
separation  between  the  plate  and  the  projectile.  The 


164 


center  of  the  plate  shows  a deceleration  (a  change  in 
slope)  at  700  microseconds  which  indicates  a second  impact. 

A second  impact  was  observed  in  plate  impact 
experiments  conducted  by  Takeda,  Sierakowski,  and  Malvern 
(1981).  They  recorded  when  the  projectile  was  in  contact 
with  the  plate.  They  observed  a first  and  second  impact 
similar  to  the  results  of  the  finite  element  analysis. 

The  deflections  of  the  center  of  the  plate  shown  in 
Figure  10.41  are  for  a plate  with  a side  length  of 
0.16  m.  The  length  of  the  side  was  chosen  to  fit  the 
period  of  the  experimentally  determined  deflections.  The 
period  of  oscillation  of  an  isotropic  clamped  plate  is 
approximately  proportional  to  the  square  of  the  length  of 
the  side.  In  the  experiment,  the  length  of  the  side 
measured  between  the  clamped  supports  was  0.14  m.  The 
computed  period  of  oscillation  of  a clamped  plate  with  a 
side  length  of  0.14  m shown  in  Figure  10.6  is  less  than  the 
period  of  the  experimental  curve.  This  difference  is 
caused  by  the  boundaries  of  the  experiment,  which  are  not 
truly  clamped.  The  effective  length  of  the  side  of  the 
plate  is  larger  than  the  distance  between  the  supports.  If 
the  plate  is  allowed  to  move  at  the  tip  of  the  clamp,  the 
effective  clamped  point  is  not  at  the  tip  of  the  clamp  but 
at  a point  inside  the  clamp.  Therefore,  a side  length 
larger  than  the  length  between  the  clamps  was  used  in  the 
computer  program. 
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The  deflection,  normal  stress,  ax,  and  average  shear 
stress,  t , along  the  x axis  in  a 15-layer  plate  for  time 

A 

from  250  to  1250  microseconds  are  shown  in  Figures  10.44, 
10.45,  and  10.46.  The  impact  is  assumed  to  be  Hertzian 
with  H equal  100  x 10^  N/m^*^  and  the  projectile  velocity 
is  19.8  m/second.  The  deflections  are  close  to  maximum  at 
time  equal  500  microseconds  for  positive  displacement  and 
1250  microseconds  for  negative  displacement  but  the 
stresses  are  not  a maximum  at  the  same  times. 

Three-Layer  Plate  with  Hertzian  Impact 
The  deflection,  normal  stress,  ax,  and  average  shear 
stress,  txz,  along  the  x axis  of  a 3-layer  plate  for  time 
from  10  to  50  microseconds  are  shown  in  Figures  10.47, 
10.48,  and  10.49  and  for  time  from  50  to  250  microseconds 
in  Figures  10.50,  10.51,  and  10.52.  The  3-layer  plate  has 
the  same  material  properties  and  dimensions  as  the  3-layer 
plate  discussed  in  a previous  section.  The  impact  is 
assumed  to  be  Hertzian,  and  the  projectile  velocity  is  37.3 
m/second.  Figures  10.47  through  10.52  for  Hertzian  impact 
can  be  compared  with  Figures  10.28  through  10.33  for 
inelastic  impact.  The  displacements  and  stresses  computed 
with  inelastic  impact  and  Hertzian  impact  are  similar.  The 
Hertzian  impact  is  a softer  impact  producing  smaller 
displacements  and  stresses  at  short  times  after  impact, 
such  as  10  microseconds.  At  longer  times,  after  impact, 
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Figure  10.44.  Deflection  of  a 15-Layer  Plate  Along  the  x 

Axis  for  Time  from  250  to  1250  microseconds 
for  Hertzian  Impact 
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Figure  10.45.  Normal  Stress,  <*  , in  a 15-Layer  Plate  Along 

the  x Axis  for  Time  from  250  to  1250 
microseconds  for  Hertzian  Impact 
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LEGEND:  TIME  e--e  o 250  500 

» » m 750  4 — i — h 1000 

z-  z z 1250 

Figure  10.46.  Average  Shear  Stress,  x„z,  in  a 15-Layer 

Plate  Along  the  x Axis  for  Time  from  250  to 
1250  microseconds  for  Hertzian  Impact 
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LEGEND:  TIME  b-b  □ 10  a 20 

*-■*  * 30  h — i — h 140 

z-g  z 50 


Figure  10.47.  Deflection  of  a 3-Layer  Plate  Along  the  x 

Axis  for  Time  from  10  to  50  microseconds  for 
Hertzian  Impact 
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LEGEND;  TIME  z z-s  IQ  a a- a 20 

■ n 30  i — ' — f*  40 

z z~g  50 


Figure  10.48.  Normal  Stress,  o in  a 3-Layer  Plate  Along 

the  x Axis  for  Time  from  10  to  50 
microseconds  for  Hertzian  Impact 
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LEGEND:  TIME  B-B-B  10  a &■  a 20 

* m » 30  h — i — h 4=0 

z z g 50 

Figure  10.49.  Average  Shear  Stress,  t z,  in  a 3-Layer 

Plate  Along  the  x Axis  for  Time  from  10  to 
50  microseconds  for  Hertzian  Impact 
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Figure  10.50.  Deflection  of  a 3-Layer  Plate  Along  the  x 

Axis  for  Time  from  50  to  250  microseconds 
for  Hertzian  Impact 
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LEGEND:  TIME  B-  Q-  □ 50  ft-a-A  300 

#— »— » 350  -i — i — h 200 

250 


Figure  10.51*  Normal  Stress,  <j  , in  a 3-Layer  Plate  Along 

the  x Axis  for  Time  from  50  to  250 
microseconds  for  Hertzian  Impact 
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LEGEND:  TIME  e-o— b 50  a— a-  a 100 

» 150  -i — i — h 200 

z z z 250 

Figure  10.52.  Average  Shear  Stress,  txz»  in  a 3-Layer 

Plate  Along  the  x Axis  for  Time  from  50  to 
250  microseconds  for  Hertzian  Impact 
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such  as  200  microseconds  both  types  of  impact  produce 
roughly  the  same  displacements  and  stresses. 

Discussion  of  Results 

The  deflections,  average  shear  stresses,  xx2  along  the 
x axis  and  xy2  along  the  y axis,  and  normal  stresses,  ax 
along  the  x axis  and  oy  along  the  y axis,  were  presented 
graphically.  Examples  were  shown  for  an  isotropic  plate, 
and  for  3-layer  and  15-layer  composite  plates  with  simply 
supported  and  clamped  boundary  conditions,  and  both 
inelastic  and  Hertzian  impact.  The  deflections  of  the  15- 
layer  plate  agreed  well  with  deflections  determined 
experimentally.  For  times  after  impact  from  10  to  50 
microseconds,  which  were  short  compared  to  the  period  of 
vibration,  it  was  found  that  the  boundary  conditions, 
simply  supported  or  clamped,  and  the  impactor  mass  had  only 
a small  effect  on  the  deflections  and  stresses.  The 
thickness  of  the  plate  had  a large  effect  on  the 
deflections  and  stresses.  The  displacements  and  stresses 
were  proportional  to  the  projectile  velocity.  The  type  of 
impact,  inelastic  or  Hertzian,  had  a small  effect  on 
stresses  and  displacements. 

For  longer  times,  100  to  500  microseconds,  the 
boundary  conditions,  simply  supported  or  clamped,  had  a 
large  effect  on  the  displacements  and  stresses.  The 
results  for  the  Hertzian  impact  showed  that  the  projectile 
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impacted  the  plate,  lost  contact,  and  then  impacted  a 
second  and  possibly  a third  time,  which  agreed  with 
experimental  observation.  The  impactor  mass  and  velocity 
have  little  effect  on  the  velocity  of  the  flexural  wave  for 
time  less  than  50  microseconds.  The  velocity  does  depend 
on  the  thickness  of  the  plate.  For  an  anisotropic  plate, 
the  velocity  propagates  faster  in  the  direction  in  which 
the  bending  stiffness  D is  the  greatest. 

The  graphs  of  normal  stresses  and  shear  stresses  show 
evidence  of  a large  flexural  wave  produced  by  impact  on  the 
plate.  The  flexural  wave  begins  near  the  point  of  impact 
and  propagates  outward.  The  flexural  wave  causes  a point 
of  maximum  normal  stress  and  shear  stress  to  begin  near  the 
point  of  impact  and  propagate  outward. 

The  motion  of  the  points  of  maximum  normal  and  shear 
stress  may  cause  delamination  cracks  to  propagate  in  the 
plate.  Experimental  studies  of  impact  of  glass-epoxy 
composite  plates  by  Takeda,  Sierakowski,  Ross,  and  Malvern 
(1932)  showed  that  impact  may  cause  a delamination  crack  to 
propagate.  They  observed  a delamination  crack  beginning 
near  the  point  of  impact  and  propagating  outward. 

The  velocities  at  which  a delamination  crack 
propagates  are  roughly  half  the  velocity  of  the  large 
flexural  wave  computed  by  finite  element  analysis.  Even 
though  the  velocities  are  different,  it  is  reasonable  to 
assume  the  flexural  wave  is  the  force  driving  the 
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delamination  crack.  Delamination  causes  the  flexural  wave 
velocity  to  decrease  because  delamination  decreases  the 
bending  stiffness,  D,  of  the  plate.  The  finite  element 
analysis  is  based  on  the  assumption  that  the  material  does 
not  fail,  and  the  computed  flexural  wave  velocities  are 
higher  than  the  observed  flexural  wave  velocity  of  a 
delaminated  plate. 

Takeda,  Sierakowski,  Ross,  and  Malvern  (1982)  also 
concluded  that  the  delamination  crack  velocities  are  little 
affected  by  impactor  velocity  3nd  mass.  The  same 
conclusions  for  flexural  waves  were  obtained  by  this  finite 
element  analysis  in  glass  epoxy. 

The  interlaminar  shear  strength  of  the  glass-epoxy 
laminates  used  in  experimental  work  by  Takeda  et  al.  (1982) 
is  not  known  for  dynamic  loading.  The  interlaminar  shear 
strength  under  static  loading  was  estimated  to  be  17.25  MPa 
by  Sierakowski,  Ross,  Malvern,  and  Cristescu  (1976,  p. 

27).  The  graphs  of  shear  stress  show  a shear  stress  near 
the  point  of  impact  greatly  exceeding  17.5  MPa  indicating 
that  plate  failure  would  occur  at  the  point  of  impact.  The 
graphs  of  shear  stress  also  show  a point  of  maximum  shear 
stress  moving  out  from  the  point  of  impact.  For  projectile 
velocities  of  19.8  m/second,  the  computed  shear  stresses 
are  less  than  17.5  MPa,  indicating  that  delamination  would 
not  occur.  See,  for  example,  Figure  10.9  where  the  maximum 
average  shear  stress  does  not  exceed  7 MPa  except  near  the 
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point  of  impact.  For  projectile  velocities  of  37.3 
m/second,  the  shear  stresses  are  generally  a little  less 
than  the  17.5  MPa,  indicating  delamination  may  occur.  See, 
for  example.  Figure  10.30  where  the  maximum  average  shear 
stress  does  not  exceed  15  MPa  except  near  the  point  of 
impact.  Experimental  data  tend  to  agree  with  these 
conclusions. 


CHAPTER  ELEVEN 
SUMMARY  AND  CONCLUSIONS 


The  finite  element  program  developed  in  this  study  is 
based  on  the  Yang,  Norris,  and  Stavsky  (1966)  plate  theory 
(YNS),  which  includes  shear  deformation  and  rotary 
inertia.  The  results  of  this  study  show  that  YNS  theory  is 
sufficient  to  calculate  deflections  which  agree  with 
deflections  determined  experimentally  in  impacted  layered 
composite  materials. 

The  variational  form  of  the  plate  equations  of  motion 
was  derived  from  Hamilton's  principle  and  was  used  to 
derive  the  finite  element  plate  equations  of  motion.  The 
differential  form  of  the  plate  equations  of  motion  was 
derived  from  the  three  dimensional  equations  of  motion  of 
an  elastic  body,  and  used  with  Galerkin's  method  to  derive 
the  finite  element  plate  equations  of  motion.  Both  methods 
produce  the  same  finite  element  equations,  but  the  algebra 
is  simpler  for  the  variational  method. 

Newmark's  method  was  used  to  integrate  the  finite 
element  equations  with  respect  to  time.  The  diagonal  mass 
matrix  was  used.  Two  types  of  impact  were  tried,  an 
inelastic  impact,  where  the  projectile  hits  the  plate  and 
sticks,  and  a Hertzian  impact,  where  the  force  between  the 
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projectile  and  the  plate  is  governed  by  the  Hertzian  impact 
law.  Both  types  of  impact  produced  deflections  at  the 
center  of  the  plate  which  agreed  well  with  deflections 
determined  experimentally  while  the  projectile  was  in 
contact  with  the  plate.  The  Hertzian  impact  method  allowed 
the  projectile  to  lose  contact  with  the  plate  and  impact  a 
second  and  possibly  a third  time.  The  computed  times  for 
first  and  second  impact  followed  the  same  trends  as 
observed  experimentally. 

The  finite  element  program  was  verified  by  solving  the 
following  three  problems  which  have  analytical  solutions: 
static  deflection  of  an  isotropic  plate  with  a point  load 
at  the  center,  cylindrical  bending  of  an  isotropic  plate, 
and  free  vibration  of  an  isotropic  plate.  The  finite 
element  solutions  converged  toward  the  analytical  solutions 
quickly  as  the  number  of  elements  was  increased. 

The  finite  element  program  was  used  to  solve  plate 
impact  problems  which  do  not  have  analytical  solutions. 

The  deflections,  normal  stress,  and  average  shear  stress 
along  the  axis  resulting  from  impact  on  the  plate  were 
shown  graphically  for  isotropic,  3-layer,  and  15-layer 
composite  plates.  The  bending  stiffness  of  the  plate  and 
the  velocity  of  the  projectile  had  a large  effect  on  the 
displacement  and  stresses.  The  effect  of  the  boundary 
conditions  on  the  displacements  and  stresses  increased  with 


time . 
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The  normal  stresses  and  the  shear  stresses  show 
evidence  of  a large  flexural  wave  which  is  produced  by 
impact  of  the  plate  and  propagates  outward.  The  flexural 
wave  may  be  the  cause  of  experimentally  observed 
delamination  crack  propagation. 
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